{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "87acfa23",
   "metadata": {},
   "source": [
    "# *eht-imaging* Polarization Tutorial\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "da8f45d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to eht-imaging! v 1.2.8 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import ehtim as eh\n",
    "from   ehtim.calibrating import self_cal as sc\n",
    "ttype='nfft'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13daf6c0",
   "metadata": {},
   "source": [
    "## Load an image "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "82370db9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading text image:  ../models/avery_sgra_eofn.txt\n"
     ]
    }
   ],
   "source": [
    "# Load the image and the telescope array\n",
    "im = eh.image.load_image('../models/avery_sgra_eofn.txt')\n",
    "eht = eh.array.load_txt('../arrays/EHT2017.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e0430b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create directory for results\n",
    "outpath = './tutorial_results/ehtim_tutorial_pol'\n",
    "if not os.path.exists(os.path.dirname(outpath)):\n",
    "    os.makedirs(os.path.dirname(outpath))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0688077b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEZCAYAAACO4n6tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAAsTAAALEwEAmpwYAACSoElEQVR4nO2dd3gc1dXGf2e1WlXLsizLvWNjeu/9C6GTQEgAE3rvvYcWSOi9GkggEHoKvfcOsY2pxg2wDe5NllVXu3u+P+6MdnY0u9pdrZo17/OMZub2Ge28c+bcc84VVcWHDx8+fPQ8BLp6AD58+PDhIzv4BO7Dhw8fPRQ+gfvw4cNHD4VP4D58+PDRQ+ETuA8fPnz0UPgE7sOHDx89FL2GwEVkpIg8LyIfiMjbIvKeiJyWw/bfFJH70ih3tIi8IyLvisgUETnUkVcqIrdYeR9Y+f/nqv9HEZlq5d0iIpKir01F5FMR+VhEXhSR/m2MbUcRec26Nx+IyCcicqGI9LPybxGRxSKyRERucdR7Q0QaRWSmiBzd1j3oDhCR40TkS+s6PxaRbdKoM05EmkVkV0faSOv/OUVEjnCki4h8KCJDO+YK0oM1jptEZLL1uzmijfJJf18icreIfCQiT4lIyJF+loic25HX4SMJVLVXbMC7wKmO812Bb3LU9gBgGbAEyGuj7GpguHW8AdAAbG6drwO8BwSt82Os8vnW+YbAYqAS8/J9FzgtST8hYB6wm3X+Z+BfKca1p1V+giNtI2AlsJ8j7R/AYx715wLHd/X/Oc3/14ZAxL5W4FhgISBt1HsSqAN2daQ9AuwB9AMWAaVW+uHAtd3gWk8G3rF+L5XWb3TjFPfF8/cF7AY8YR0/ARxpHfcFJgOhrr7W3rj1Ggkc2BpDjgCo6nvA4zlq+yDgPMyPedc2yl6lqj9bY/gO+A74lZW3EEOCEev8f0AZhhwAjgNeUdXlqhoDHsI8oF7YG4iq6rvW+d+A34nIAHdBEQkA9wF/VdUZdrqqfgPcBTS3cU09DesB1Y5r/RQYDAxKVkFEtgJqMS9qJzYHPlTVVcB8YLyIFAJnA9fneNzZ4CTgH6oaU9XlwEvACUnKpvp9bQ58aB1/CGxpHV8K3Kaq4Q4ZvY+U6E0EPg+4UERK7ARVbXnARGSs9cn7qYg8LiL/FZG5InKiiFxhqQ7uEpEnLFXBPxxt7w48BbwKHJJqEKp6myupEIsUVLVeVedY4wkCxwPPq+pSq+xWwAxH3enABiJS5NFVQlnrpVGPeRDd2BwYDbzhMd4rVfX1VNfkhojsKiIzLFXMeyKyQESmZ9KG1c7WlppjrohcYKk6/icio0Rkkoh8LSKPZNou8DGQJyLbW+f7Y74gliatYb5g/uyRHgFsNUPAOj8TeFRVa7IYGwDWby4sItNEZF0RWV9EvrPua9IXjauNAmBjWv9mtvSukfL31eo6RWQEsAPmy8RHFyDY1QPoRJwBPAMcICLPYqSS9x35TwDPqep1lt5yOkayeABARMZgpOttgHyrPURkILBUVcMi8iRwj4ic6pCik8Jqsz/wH1f6ERjp7UeMdG9jIEalYqMa81BVAj+7mneXtctXeQxlrLVf0NaYLfxaRN5zpTlJJQJcrapPWGTzDXCOV0MisjXwLDBWVRudear6PxE5G/Ni+VBVbxKR54B/Y75aGoGfRWRbVf3M6uupFOO+XlVfU9WFInIQ8B8RqQFiwJ6qGk0yxr2B71T1F48phw+s+/EZRpW2HDgY2C7FONqEql4tIpsCU1R1pjWO14H7VHWxiFwM7JWk+mJVPZS4KsT9m/H6DUDq39eHwOXAvZh7/wjwV+AKVfXjcXQReg2Bq+rblsQwEfgj8K6IPKiqJ4nISIyK5fdW2QUi8r5HM2+par11/BdrfxDmxQDwIvB3zA88pdRqTQ7dAZyoqmtcY/0n8E8RORb4TEQ2dUhzXg9LsonMTMom5InIFsAtmIf3NVU931HuTVU93HU9cx3j/8iR9TeM7j3Z/VgNfE9qNc0aVf3MOv4WoxpabfU7CxgDfKaqi2lbhYWITMAQ/e6q+rWInADcLyK7u0ncUi9dBPwuSXNXAjcBRwFHWmVvBCpF5AagFCMsvNDWuDzwKOZFfq2I5AGjVXU2tHw9pquicf8OUv0GPH8zqvqFiLwlIv8BvgJ+Afqp6jsicglGqv8JuERV1zaVW7dFb1KhoKp1qvo3Vd0NMylzvCUFD7aKLHcUX+nRhFuiBfgN8GdLIn0Vo6Y4OI3h/BV4V1VfTDHehzAPlN3eUqDcUaTcynfrZb3K2uW91AQ/WvsWiwlVnaqquwJfYkg8Y4jIicC6wAXJyqjqTFVtRZwuOF9wEY/zEJnhGIxk+7V1/jdgC4wqzI3DgNdV1ev3gKpWq+oJqvo7zBzG5qr6DObl9yqG2O8WkfIMxwjwMuZFsA1mkrmViqsNLMd8XTj7Lie5qijl70tV71HVg1T1z5iXx8UisgdmUvcgzBzQ0RmO0Uc70GskcBG5T1VPsc9V9X0RWYH50S2ykgdgJqLAqDbmttHmIGCOqp7uSDsEuE9ETk4miYjImRjLklut83GqOltEtgNqrclDG3WArbefjCFEG+tjPu0bPLqZDDhNFIcDxcBUj7JfWNe6F2Yys90QkbEYSXRvVa3LRZtp9JmWCgVD+C3/G1VVEYlgJozd2AnYUET2tM4HAbeLyDxV/a2r7HXAn6zjTTET1mtE5BdgPGZSOm2oarOIPIV5CZRgJkYBSEeFoqpNIvIN5jdj970+5rfhhbR+XyKyL/Czqn4rIvsQ/01NBjZL9/p85AC5Mmfp7hswG9jacb4LxmSqyDr/HPiTdTwUI4Ff5Sj/D+e5lXY6hqCcaWVAE7BPknEcitF598F8Xtuf2GCkl/uxzNmA7a22NrHON8S8bPpjvp7exmFGiFHl7GIdF2Ambu3zK0ltRrgX5rN4Q0faUMyD/5DrPqQ0I7TG9hHGqsXOn5Sk3wnAWyQxv8SoROY6zq+y75d1/h5wdIa/hX0x+t3B1vk+mC8n27zzj8CtSerOxWFG6EjfDjPhbJ+/CvzW+j8ssPvK4ne7tfVb/GeW9U+2fidi/W4WO35Ptm67Ip3fl1UmD0PUQ63zQ4BnrePbgQuyGae/Zbd1+QA67UKN6dT7GNvWD6wHf1tH/liLdD4FHrbI8Eor71zrhz8XuNFKO8RK+8RJPhizuzqr7G9cYyjBELK6tn9Y+SMwplsfWWOcAhzsauOPGIlnMuYzXRx5XwMHOc43s67nI4x+vn8b92gn4E3roZ6MkcwvBsqs/Fusa14C3OKo9wZmQnEm5iV0hHVdz2MmHP+Ng4RdfW5rEVyhR976GBVOI+bF9hvrvi4GTgGuwBDxDOD/Mvw9nANMc1zrvo6883GQsZVmm6E2WmO6wpX/PrCe43wTjFAwGTg3xTg+AX7Xxli/B36b5e9eMDr6ydbv5ghH3nCMemRYOr8vx3N0reM8iPnq+Rh4xf6t+FvnbLak1+shIhXq0HOKyMvAS6qaE5WCDx9uWPbiPwMbqZmA9eEjI/SqScw2cLtlnWDri7fHfEL68NFRuBK4zidvH9mi10xipoHXMKZ7dRi99MmqOquLx+Rj7cb1aplD+vCRDbpMhSIiFwCjMKZO4zBuvEXEHVjGAZeq6hKPuodj9LtR4AdVvb+Thu3Dhw8f3QZdQuCWudd0oFJVYyLyPGbScCfgHVV9RkT2x0zgHeGqOwwTz2EzVVURmQwcppaDgw8fPnz0FnSVCqUeCGNM7qoxKovvMNL3X60yH2Pcdd3YE5iq8TfPp5jATa0I3HIkOdE63SJHY/fhw0cSqGoqL09P7LXXXrp8+fK2CwJTp059XVWT2b/3OnQJgatqjaVCeVpEFmHsj+dgYjTYXnY1QD8RCWpiXBFnGbucZ2wHNXFM7FgmaX1q9LRZ3Z42XjdiXT2ALNATx9ydsXz5cqZMSc/HSSQvK6/gtRVd8vxbQXouwNjeHo3Rg1+BceXtYxUrA1Zp66BQzjJ2uVRR5NJGTyDDgGvr6eiJ19NTxtlzoJjXYjqbDye66rc4FFjpIOdFmLCqLxOP4raDdY6IBKxAVGCCRG0h8bBw22G83tqF7vxQ9jSCaw96yrV29/H1PPgEng26Sgf+GrCPmGW5qjEuvGdjvBRvEJHxGM9IOwLexsA/MQ4Pv4jIzcBtIhIF/tbeCczu+jB213F1FtzX390eX3t83W1cPRP+XcwGvcYT00sH3h0JsjuOqTuiuz3u3W08XYVsJjG33HJznTLlg7TKivSZqqrJFqTodei1jjzdjSi723i6O7qb9Bug+4yl58HWgfvIFL2SwLsTWXansfREOO9fV1NAd3up9Cz4dy0b9DoC7y6E2V3GsTahuxCoL41nCl8Czxa9isC7A2l2hzEkQ6Zj666PXHcgcp/EM0Vu7pYV4XESJs58EBOOw3MlIxE5D7PiEsBTqnqTI28UJqx0EBOy4xhV/cnK2x+zJm4Qs2rRP1T1Tkfd1zBWdTZqVXW/XFyfG72KwLsS3YG4cz2GdNrrahLtyjH4JJ4JcnanrsIYZ2xrWbN9JiLruWMqichemNjmm1pJX4rIdFV92Tp/EnhAVR8WkWOApzEx4QFuA/6gqtOssCDfWSs0PW/lL7b8Wzoc3YFX1np0xU12O8h01T/aaxydPZauvn4fbUExS5umsyWHtQD18ZiFxbGiiU4DDvcofhLwhKo2qmoj8Dhm9SJEZBMMsT9ulX0c2Mha6BvgHlWdZvWxGLNIzB4ZXnRO4P++OhCdSRzdgawzQVeMt6vuTU/4f3QtcuaJOQazHNwMR9p0wMvscKsU5bYCflLVMIC1n2Pnq+ptrrYKSVxYvEREHheRD0XkOeuF0CHwVSgdhM54aNcmYuhMp52uUK346pS2kPbdqRSRKY7zB6yYRwADrb0zxno1Zmk+NwZ6lKtKkufOb4GIlGEI/zRH8g8YB8M5lr78IxEZr6qL3PXbC5/AcwyfuHODzjAP7Gwi90k8FdK+M8vTcORxO+0lcy5K5cXolefVzg3A1ao6r6Wi6sWO4xdFZDpwpFU2p+gNXNBp6Mib2ZPUI7lGR197Z97T3vj/axs5U6HYQe3KHWnleAe7W+pRblmSPM92rHDVzap6Txvjmg+MbKNMVvB/TzlCR5KL/0+Ko6PIvLN18T6cyM0kJkZ1sRJY15G2PjDZo+zkFOUmA6NFJARg7cc62xGRA4HdgLOs83HWvkpEjnb1NRBY2Nbgs4H/W2onOurB72riTmY90h2sStxjzHWbnQH/wXOj/RK4qsaAB4FjoYVUNwUeF5H1RORtEcmzik8CJopIoWU7fpiVhqp+CXwFTLTKTgS+U9WpVrs7A2cCp2MmLEuBy6yyxcB5IlJild0KoyN/Opu70hZ8HXg70JM/6Tub+HrSpGRn6cZ9nbgTObUDnyQin2H4baKqLrYccyYA+UBUVV8TkQ0wK38BPOSwAQc4FHhIRI7HOPIc7Mh7EhiCWcfAxvvWfjHwHPCmiESAEHBgRy352KuiEeaStDpK6u4odEeJryPJK1dtdzTBrm0Enl00wg10ypT0BFSRjfxohA74EngW6Amf7d2RsN3oSNPBXEnRHS2N+1K4Df8uZAOfwDNEdyfvnkDcydARhJ4rguxIovVJ3A9mlS18As8A3ZVsezJpp0KubMFzKY37JN5RaNPCxIcHfAJPE91Rf762ErcXckHCuWqjdxNtR8CXwLOFT+BpoLuRd0/UmefaQqQ9bbaXyDuKxHv3y6H3Xnl70GUELiLrYuwrG4BdMOY/S4HLMYFjRgHnqWqtR90LgDKgH/CGqr7QUePsTuTdncaSi/66eoKxPYTpk3gu4Uvg2aJLCNwypr8V2F9VYyLyKEYJ9k/gClX9n4icAVyEIXRn3W2A3VR1HxEJAt+LyPuq6g4+0250F8LsrlJ7e5GrScv2ELlP4t0Fve+Kc4Gueq63wgSGOUNELgH2x0T72o24u+rHwL4edfcDPgVQ1QjwPUaCzynWBvLuak/JTNFe787OrmfX9ZEL5CQWSq9DV6lQRgLbYbykVovIY5g4vg0a9yyqwSN8o5X2veM8WTk72MyJmQ6uO5B3V0vs3QHZ6ruzlcizlXw7QmLuXVK4HQvFR6boKgKvAWY41B4fATsBRSIiFomXkTyKWB/HebJyWHGCHwDjiZmjsaeFriDg7q5q6WwLkmzrdBcS7z3wdeDZoquEtc+B/o7AMiOB7zBLE21lpe0AvAwmGpiIDLXSX8ZI74hIPrAe8EGuBtbV0mtXqA46S9WSi/6yqZtpP91FndLVv8XOha9CyQZdIoGr6koRuQi4XUSWAQOAq4EngCtEZA9gBHCuVWVP4BRgH1X9TETeFZFrMVYo56lqdS7G1dX22Z1BNN2NFNozkZmJhJ2pNN7Zahgf/l3LBn4wKwtdSd4dTdydZTfeEY9gJm2mWzbTcWZzXbm8Fz2J2rILZjVOp0y5I62yIvv6wawc8B15WHvJu7NVAenWy0bKTqdeulJzNtJ4V0ria79U709iZgufwHOEjibv7qjzzRbZqk0yIeh01SodSeI+MoF/d7NBryfwrrLc6AhCznW5zkKmhJ4Okeea7DMtm035zmqre2LtvrqOQq8m8LWFvDuSuHNJ9plK2m3VyRWRdyQ5rv3Emwv4ZoTZolcTeHvRU8i7u0jmydpPh4BTlcsFSXeExJ5rrN0vg7X3yjoSvZbAu0KNkEvy7mxy70jkSuJub75dJtckvnYTby7gS+DZolcSeFeoTnIpBbeX3DvDkiUZ0tVvpyqbCyL3Sby7wbdCyQa9ksDbi46SWjuavLuDRJ6J/XhbZN4eou4KEs8F1s4XgS+BZ4teR+DtJaiO0nu3t0x7iLsjzQ4z8ZJMVT4VWbcnrzMnN9dO8s0V/DuTDXodgXc2uit5d5Y0nunEZXuk7mzIur2k6kvhuYAvgWeLXkXgnS1950odkUuCzrU0ni3Ssf1ORea5JOv2qFrSLZNN2d4F/65kg15F4O1BV01aZkPE3Z3U2yJrrzLJCDub9I4g8c5GdxxT9vAl8GzhE3gXItfk3dHpuUI6E5kdScwdQeK+FN5e+FYo2cAn8DTQEdJ3Z5B3Lsq2lZctkhGzO6+t9LbS7PS1VRJfO+BL4NnCJ/AcozuTd3tJPpsy6UxWusu1l7TTJexcT2z6Unh74N+NbOATeBvIRPrsSKuN9pB3Z5B5MrRV135sc03a7SH2bNLTzc8l1p6XgC+BZwufwFOgI1QHHT352BZ5Z0v46eSli7ak7VRp7SXtziBxH9nAv6PZoLOsxtZ6tFfl0F7yDiRJS5YfcKW5z91pycpmumXSbjrjbOu60zlPltbR6Mgvtp4Fe0GHdLbUEJFCEfmHiHwmIlOs5RmTlT1PRKZa2wWuvFHW0o0fish7IjLakbe/iLwhIu+IyBcicqarboWIvCAiH1nj2DztW5EhfAk8CXL9YOSyvfYSlpukU5VNlpYq3Qte0rUzz0sn7paovc7TyUvnPJO0VOlt5flIhpzdsaswS0VuKyLjgc9EZD1VXeIsJCJ7AScAm1pJX4rIdFV92Tp/EnhAVR8WkWOAp4GtrbzbgD+o6jQRGQR8JyLzVPV5K/9eYJqqXikiuwHPi8g6qtqUq4u00aUvcBEpEpGvReRm67xQRO4WkUtE5CHrH+BVb3cRuVdErhKRK3M9rkxvSntvYkepSJJJ1MmOcyFhJ9uC1tYeKT3V2NO9flKcZ5LWFtYOybizYOvA09mSQ0QCwPHA3wFUdRYwDTjco/hJwBOq2qiqjcDjwMlWO5tgiP1xq+zjwEYisoV1fo+qTrP6WAy8C+xh1a0A/uAYw7tAGNgvzZuREbr6d/YXzA22cTYwX1Wvw7zl/u6uICLFwCTgHFW9CthYRH7V8UP1Rjo3MFWZjiRv53Fb5NcZJG4TeTppmZJ3JqROivNMkIuHp6sfwO6D9hM4MAboD8xwpE0HvBZB3ipFua2An1Q1DGDt59j5qnqbq61CYJl1vDnQpKrz0xhDu9Flvx8ROQL4GPjJkbwv8CmAqn4DbCIiZa6q2wHzHJ8jH1v1uiW6A3m7j9Mh8Wyl6aDHlqqMF4mnS+jpXE+y+5ApqXdnou3OY0sfaRN4paXbtrcTHY0MtParHWnVQJVHhwNTlHPnJW3H4qetgIczrZsLdIkOXETWB9ZT1UtFZGNHVhWwxnFeY6XVpFHGq58TgRO98pIhk4ehIx6cTMk7WV4y4s42ra1xpIJXOS/dtfvYrc92pnm151U/WZteY2xLR55OvXTzsim39iIjM8LlqtqWNKuuc0mzXFt5Xu3cAFytqvOyqNtudNUk5oFAo4hcDOwIhETkbGAp0MdRrsxKcyKdMgCo6gPAAwAikuqfBeSevLORvtsql6mknSytrX1bacnOM0GsjWPnPh0ST1U2GaHjUYY0z33kCKoQyYkrvc0D5UmO3WXLHeflxNUg7jzPdiwBsVlV73G129ej7nepBp4tuoTAVfWv9rGIFAKlqnq7dbwd8KGIbAR8pao1VrnRqvoTRsUyUkQKLDXKDphZ326FXKhOckHe6eyzIfV0xtsW2iJte+8+Dngcp9OXF6G70VZ7Xvm5kMLbix7/conlZPQ/ACuBdYmT7frAKx5lJ1vlcJSb7MgbLSIhVQ2LSAgY68hHRA4EdgMOs87Hqeps4AugUESGq+rPjrb/0f7La40uVZ+JyEHAzsC2IjIRuANDzpcB5wHHWeUGAB+JSKGq1gOnAHeKyF+Ar1X17c4cd2epTtLJT5e8nVvQsQ8mSfPaQo69+zjZVujavPJCHnuvvoKuY6/r8Nq87kmyNDfWDv1yN4cC0Vh6W6pmVGPAg8CxYEgVy5pERNYTkbdFJM8qPgmYaFm+FWKIeJLVzpfAV8BEq+xE4DtVnWq1uzNwJnA6UCIipcBlVt0VwL8cY9gF87O1zRNzii61A1fV/wD/cSWf5lFuGTDUcf4m8GYux5LLBzVTQmgrLRfStpvM0jn3Ok61t5V89hOS7Hrcj2HU2juNyXAcZ7Olg7bULaR5niwtE/R4Cbpd0FxJ4ABXAZNE5DMMv01U1cUiMgqYAOQDUVV9TUQ2wBhBADzksAEHOBR4SESOx/xED3bkPQkMAZY70t53HJ8K/ENEPsI8Dr+1TBVzDlFtUzW8VkBENNXbKluddCZlcqk6yVQ1kkoyTTcfDEnnufLsdFx5XtcitJ7hcZJulEQij1qbfR4hkajd56kI3X2cap/q2Os8WVqq9GzLdVT9XEBVM56s23LzQTrlAy9T7daQPrdMTWMSs9fA/0Kk86TvdMq3l7y91CX2VpuXx6hdd6UmEEiq/rDVGqsqKijccksagGKgCCgFSqxz+7gUWFJSgmy7LTWY2eVy19bPsS0qKiKy9dbUONIqrM1ZZw1QsMMOrCwtpcxq19l/sTXOACBDh6IjRpBPaxVLKnPFVPcw2b1Pdp4sLVV6rtFjH2a1JPB0Nh8J6LH/865Ae25WJg+8V34uJO7Cvn159d13iYq0IvBCoMBxPGGbbbhx0iTAkLcXcdukOnD99bnv6aeBRBKuwBB0f2CAtVWstx6T/vtf1JFWaZWx9/0wUvgDzzxD1bhx9LH6KnHs7ZdKMXDEGWdw6nnnAd568lQ26ckI233fk+XnEr32YVSMFUo6m48E9PpYKF0pfaeq317J20sSz+vTh4aGBgLRaIv0modRCuZZWwvZFRTQ1NRECcYmSqxyTinWrtNYWcny5csZCAwOQDAIAcfe3iQAC0cM45dffmF0HxheAOoQriIRs2+OwLwINDU1URgKMRCjKglbWxNGrRLGvHR+mDWLA3//ewRD6M7QR241S6YUkEo/nokuvHfruNtCTnXgvQq9nsDTRUdL39nqzdMhbntfWFpKbW0tIQwZ28Tt3NvSqxYW0tTURBlGmhZXuaBF0MEgrB5gCHxoOQwpNWl5Qch3EbkIzB42nJ9//pkJQ2B4USJ5N0cgaglaFTUQDocpKCqgKgThMDRiyLsZcxzGfC0smD2bcePGEcNI5c0kkniAOIlD8glPr0lN997r/+BTTzthW6H4yBgZE7iIDMQ8JwvtWAE9Fd1J+m6rrUxUJclMA4v69GHNmjUUYNQQtrohj7j6xCbwuoICGhsbqQD6B6EgFCfsUMixz4fF/Q2BrzcCxpRDfiixvAQgLwDVYRg81EjgW46DQcVxqTsSMSRt7wcuNBJ4QXGIIUMsAm+Mbw0xQ+KNwNzZsxk9ejTk59OnubmVlJ6MzN3SuJOMsyFzH9nCl8CzRVoEbkX5+jPGLtt+BkpE5D3gLFfgFh8OZCp9t6U6yYS83VuxJYHbk4C2JF6Ay0Y7AM2WCqWyHAaXGrK2N5ugCwuhUaCif3+WL1/OZmNhdH+TFwolSuGBAExbDMOHD+fLaV9w8oZGqndK3TaBNzbCoGYjgYeKChg2DBoaoaE+TuD11nF9I+QtXkx9fT1Vo0dTNmsWYeISeh5GIg+77pNN3qnMD9MlcTehdyXB99iXi0/gWSFdCfw64DPgr057RhHZFLhGRM5S1ercD697oC3pOlP1RyZ1vSbYUpG3c7LO6QwTBEosCdyeiLQtTgqA4oAhZHtbaRH40IEwrF88PRSC4uL4/ud6qKysZOHChWyxST7FpYUQyodQgUMJLhAIsLC2jmHDhvHxOy9RuP1wiMUoiDRbzB2FcBOEwzQ1wtBaI4HnF4QYNiyRtOvqoa7WHK+phZFrYPbs2QwaN46KWbNoBOoxJB6y9kHi0riTzJNJ4zZSSd5t6cO9/qe5KLPWQdW8xX1kjDYJ3JK+7/GSslX1SyseQH9MxK0eg1yqPHLVrxdZu9Mzkbpt4i7ESNq2BF6GmZgsAoqDUFQIpaWGoIuLobAIZlkEPnwIjBkczysohJJikELD4D8uzKOyspIZ07+hePeN4sQdDEIwzxC4hV++nMeGw4cTaVoOo0ZZYndTAnkTbqagsYHyFfU0NTURLChgaH9osgjc3mprzb60FsaHDIEPGT+ekpdfbvmaaATqMFK4fW6Td8R1nC6Z2/8Hn4hzDF8CzwptErjlnjofQETWwZjoNgNHAs9aUbgWduQguxK5lL7TVZ0ky29rsjKZy7staRdbEng5UBEwhFxamrgvK4O8AiiwCHz94TBmBOQVWyJ3cQkUF1n7YlbU11JZWUm4cTWMGWMI3J65dCISYUH9QoYOHUoeq6F8ENjSd7jZEHhjk9mH8gmVCeFwmEB+IaXlQYoaIxRYL5H6+sT9hoUwfdYs1hk3jsX9oLIJauuNFJ4PNGDI2/mysydv3CqViFXG7Sjk/F+kq0rxkSbU14Fni0wnMS8HrsSEUFxqHR+b60F1NLpK+k4HyfTeyaRwL6nbrde29d22BF5VCJVlhrTLrH15efw8mh8n8MrxA8gb1ddklvaJk3dhIRQXsXzWXEZXVhLVBlhnfCJxt5iYNEMkSnVzH/Lz8yktaTLthK2yMSsaXTAPYvkQiVCQH6C2qYlAMASFReQFmigJhMm3Jk/tLwJ7zG/+NJvtttuOX8qgSqCwxkjphWGoxZB5kLg6xalW8dKP24Ruk7ktpbdcGm1Panql+ySfBD6BZ4VMCXwasACYoKpHi8iFHTCmboNspe90yrYlfaeqky5521uRlff999+zaNEi1i+DykpDfOXl8X2fMigtD7JUSvnqi6/Iz8/noN+OgpH94gReWBAn8FCI5U0Leemll6hvXAUFlY7RNkNeBIgaiTzWwJLljdx+++0Mqcgz0nZ9HdTWWcrtBqhvMPvGJkKrG3jzzTeZM2s2DaVrEnw5bFvxWMw0HQrBku8n8/bbQ2kMQkVfM4EaCkGoFoK1lh08yZ15vMypnOoUt6QN6ZN2pmhvGz3uJaH4BJ4lMiXwjYE7gTdEpAizhJEPF7KV8N0E7TxOpvfOA2TjjQmKEPnqK0oxhF1I3GOxBCjJh9Xvvkj9/+3LR/0Gs+OQRZSXQ/8KQ95FlSWGySsq6F/SB/lwBiUDN+CT1aUcNGE9SwFuO7DHLcl32SXKrbd+xsqVTZjQ7DHidBiDqKUWqa/nuN9W8PQ7H7FmaTPn7LgAamqguhpqaojWNlJTE5+snPwjVA+chy75kS+/TDQTdtuXv/gzDFh/E2a9+CAXjTcamVDI6OoLC43evrAW8hvhi5EjGbzeejS89hql9v0jkcgDGLXLurvswoJffmHNDz8kSOJeRJ6JFO7DDfXtwLNEplxzPWYdueuBbclxRMDOQFeoT9KRvr0mMNPRe8eAo084gb1///sEy5JC4uTdtxD6lUNeMfzl5tuoGLMOgwfB0CEweFiAomEVMGwojBgOI4aTN3oE/YetzylnXEi4Tz/oNwwKhgCDiDu8lwN9yMsr4be/PYQNN9yaRHmgGdSSsmtqYOVKIitWcewJZzNw6Prw448wZw5r5izl51mNzJoFs2bBjBnw/QyYvwBOPf8SBq+7GXPmwA/W9tOPMHcuzJ8Pv/wCCxfC4lVw+V9vonDAyBZva9tKxv66KC+HfsUwbtNNuehKsw52iWNzh7HNB0449VR+ve++5CW5/+n8f73QnVV4XQLbCiWdzUcCMpLArVWeZ1mn71pBzddK5HLyMhdtJFOb1APjxo3jmSefpA/xuCWlGHm4tBgqKqBvXwivgaFDhxKsW8CIzSC/si9UDTAFKgcYlquogLIyJPgzjY2NFBYNxhC37btpj8RM8xVYDj/BYD7RaCN5eWGgDprjxM3KFbB0KaFFvxAOhxHymPfFSqqrWwTwFquShkZjcfJLI4gI0WZl5o+JUqytCsmzHIwWxCAQCNBYH2P+aigqijsbFRcnOhY11qyiX79+BDBxVZKpVKLAvHnzGDlyZIuEnux/409i5gC+CiUrZETg1urv12HEMMFwxLMdMK61FsnI2Z2XjvRtT8Kts846LJwzh3LikncRhrxbpM9+QF45xcXF9I0tJH9IJVQNhKoqqOhnlOLl5VBUDvRBZLHxhCzoR9yZ3q00iFBQEDAxSwoLaWqqpbi4GZrWQPUqQ97Ll8PSZbBwAfnLVtIYDkNeAT/OipP3mtq4k07Y8rBcjUXgqiyzerTjhoMh8UDEqNlXYhF4Y4yFC63JzT5GhZIfMl6ghYXWnGptNeXl5RCAPrHW4W9tKDB/3jx22XVXophXl3si03knnMe+dUqG8K1QskamOvCJwJ6YhYJvxaya02OQq0/XTCXnTNtOReTuSUvNz2fUqFEsmzOHwRipsgwosyTv8nLD0f37Q6hwKCtWrKCqqsDYYVdWwqDBhsD7VGCsw8uAYgKBoEXM5jxuWGdH7I4ADRQVKY2NjZbVyiqKNWwRt5G6WbgQli5h1cJGli6GwnCYQDDE0qVG6q6tjXtX2nFOwtZeRIip0ujqGRKJsxZD4GFV5q2BojVQssyad7U224s0Vmck8FgBlEchL2za8QpivWj+fEaOHImS+KA4vTizJWqf1F3wCTwrZErgM1V1lYgEVbVZRPp1yKg6AJlEme8oHWW6um/7PNnEpU3geUDlyJHU19fTvGRJPORrYVz3W1lpCHzAAAjmDWPBggUMHTMAhg2Hyv5GCs+vwEjZcY1wIJBPY2MjZWXFVm/umH6GbgsKYjQ2NlNYWEjjyp8hFDYS9+JFsHgJ/PIzK5bDosWwphoC4TB5+aEWlUk4bFmUYMU0sq5fMQSOteCI87XhJHKskQQCAZpjMZZipOWSGBTVQL/auIlkaSnE1qwiFAqRX1ZMn2g9gVqMyO+BJZYKJUr8QXG+OLykcTs91cRmOmgvwfe4F4RP4FkhUwLfRUSmYhbt/BswrgPGtNYgW+nbPVGWTAIvAPqtsw5z5syhEmvSMhhXZVdWQv9KGDQI+laFyG+2CHzCYBg2zGTSn7jkbU9/5hEM5lsqlALixN1sbU3W1khhYTNNTU306dOHxsUroDgcV5ssX8bqaqMmaagHrBCxwfwQqnGHzVgsHlYWq+U84gRuqzmaif9g3Sv3iAjRWKwlwJUdsbA5BvUroaQG+pZBuLmWSCRCXt9yShrryQtY/YZbRyucP28eVVVVhIqKCDY0tORBojrFST2+NUoWUN8KJVtkSuCHYJ6dz4DjgRuz6VRExgJ/wazgPAxYoapXi0gFxsLlR8zL4VJVXeJR/3BgM2ssP6jq/dmMIxtkO3mZifTtLJdKfZIPDBs3jjlz5lAFlAYMSdl67/6VMLAK+g4qpKGiir51g1i4cCE7bDgK+g/GkHcFccnbTFJGoxAMBi0CzydO3lGc5A31FBY2xVUoK6qNLmT5Mli5kqbqRtZYKpJYDIJiAlTlh0IEQya8bMJzG4aAdR7ESNWqSr7jfjjtspud9zAQIBaLteTXO0ZZC5RGoG6labi6uppgWT/yYwsJBKx1AmrBjq3Z8mJYvZqamhoGjhhBZOZMIE7czqlc5/8vl1J4r4Gqv1hDlkibwEXkG+Brx/Zf4Cjgr1n0WwE8parPW21PF5GXgROAt1T1GRHZH7gZOMI1jmHA+cBmqqoiMllE3lHV2VmMoxU6Wk+eqnwq6dupPrFje0SAsZYEPoS4V6KtOhlQCQMGBWDQIBZSwdChQ40EfvAYzFo45dZmS96m13C4scUTs7DQXiLBGZzV9l9sICiRlknMxlVrQBuNQ064iWbrmbRttosLYGU4TEFBAZIPIWn95SyWG2SAuAQeIpG8Ia5qiWFNrwYCqNWYTcBgSNyp+IlFYNWqVeSV9KNpsRlbYaHVZo3Jt9fgLMVYogwcMYJlFoG7VSmpyLst+OTugPp3IhtkwjO7AA9ifBwOBb4F9smmU1WdbJO3Yxx1wL7Ap1bax9a5G3sCUzW+GvOnwN7ZjKMr4KXrTlYumQRegJGVbQuUeXPmMLAoTt62+mTgIMwk5aDBLKAPQ4cOZcWKpfTpMwKowrxH7YXRbCvoAE1Nlk67sZGCAlvWjZCoQmlEY41QX0csZpVfbZmSRJohFiMvQIvreygEJYVWiFhLArdd4osckQ4LQqZOEEPgotrKdNK5GIVNqIFAAFRbEar92qnHWLaEMQQeKutHdbXRwdvmhiXF5ivGtguvwBD4kJEjE5Zpc75QvV62ydBRcytrBfw1MbNC2hK4qq4E3rM2RGQccFl7B2DZkr+uqjNEpAoTLAsw695aE6ZO4ctZxi5XlaTtEzEWMzlBLtQnXmleZJ1MdeIkcQGe/uc/mT15MusXJxJ4vwoMi1uzmGtWRHjxxf+yYMFPwO4YnXeRtdlaZvNODIeb+PDDD4jFlEMOOQ5DgxHXvhlpbkabwnz5xf9YsngJp27RHH/IAoF4UEKLxItC8Morr6CRCMfkm7jjzmczP2SqxmJGlTJp0iR+mj2b4SnubbN1H6699lqqV6xgoCPPlrwhLr0rcM899zBn9mzWqTV9lZaaMeZbL5RwvSH6cuD+++9n5aJFCffebtd57CVNt1eN0mskdN+MMGukLRSIyHjnuaWy2Lg9nYvIbsBuwDlW0lKMSAiGYVa5yNtdxi631Kt9VX1AVbdU1S3TsULpbPVJWy8EL8sTW/osBCYAI556iq0X/cBrG+/Cv2MV9LVMukurSoz0PcRs+0zchhtu2IxotJADD3yQxYvDxGXNfJyjqarqy/DhxWy77bbce+/zxHXg9lo3lqFfuImG6lokr4TaugYGl1svgkAeBPNbiNveFjVB//79aaxfQ2WRIUw7Lxi0JG9rK8izJPBg0PMFZq/HmY/5fmhqaiJUVtaS53SPh7gKpQ7o27cvkcZGlgHL62lx4c8LQFExVIfgGxEGAn1feokhU6e2BARzjqGtF29XozuMIW34EnhWyOR/fL+IzBeRT0XkfhF5BPhWRIqz6VhE9sWoQ84CBonIdsDLwHZWkR2sc0QkICIjrPTXgS1ExObk7YBXsxlDZ8Nr0tJLEnceuyVxO1hVyNoXh+Cn/DweefpfyOCx9O8HfSsCcem7cgBUDkCkgpkzaznhhBMYPXpdBgwYBC3aZbfiIUZeXiFnnXUWRUW2C4utGVag2VpGJ0osHGGXXXZhr732IiCO12QggATjUritB//DH/7AZptvYdbNDCTmOTcROOaYYxg7blwLGTv3NlHb56eddhpDhg4lgDeB22gGzj77bKoGDaIB8ym32oq/EonAEoXIMSfx+7vvZgXm+8S5Vmgy8nb/L73+v+mc90rYk5j+qvQZI+3fj6rupqojMJYoLwFzML/vL0VkRiadisgWwNOYeCrvAs8D6wKXAr8WkcuA32EmK8FI+i9b4/gFM7l5m4jcAvwtVxOYqZBL9YlXmVQqFbf1iT3lWIK1EMNWO9LU1MSI6ilUVmKpTQbEt7xKoJIZM1YzYcIEamtXkZdn673dI7QmAjWP1atXU15eTJy8bV14tOWBikWiLRYggQAkRJkKBFpO84JWtNhYDERaEbsEWqq0rGBvw9Z324TtJGjBvFKCwSCxSCThntmvE7cUXlpaSkNtLY2YCZ1qoNpyJqoCZr/5CkcccQTV22yTEJbXi8idZO71v3P+T30kgb2ocTqbjwS0qQO3VuTpp6orAKyVeeYDLzryh2TSqapOxXz5euEEj/JfAhs5zh8DHsukz85Ce9Qn6erAW0g8BLUh2H6/A3nuuef43XpKUXnIKMHt+CZ9KjDa3DJmzlzMYYeNIy/P1uDaiLXaiwSprq6mb1+nF6Zji0UhFiUWjcUJPEEClxYCDwSMtJ0XgKgqIoEWwhZHftSR5qXyckrhUeJk7iRwG/YMt03wtj65GSgpKaGpro4IcR+ePCBUb/rfbtnP/OXqq7lu0iTO3XJL1o1GacR88di26Mn+Vz7FZANfB54t2uQba0Wey0RkR3eeiAwG7sP8rns0OlpCSqY+8foMdxO3U31iS4OFGMuJ76JwwAEH8L83/ssWo4CK/i1qEyOOD8DYU5SxalWUSCTCgAHlrl4hTj3GRSYQCFFdXU15eRFx4rYM9Gyvm1iMWEzjBJ4XMPpvhyidQNJ5oKqICHkO4gacK68BiRK4TdTJzlsIPBpNUJ04y7SUDQQoKSmhsba2ZUq2CWMrXg1U10JpM8yYdDt5eXlsfsYZLeRtS+Ft6cHd8NUmaUBj6W0+EpDub+lS4CgRWSgi34jINBGZB/wHs15mK2eb3oBcqE+SlXcSgm0uZ5OHLX0XF0PD+M0pLS1l4JIPyC8vgfK+cQlcyombCfZBtZCZM2ey7rrr4C3jxh+QUKjQIvASV140YUKplQoFWqRvAnktaTanq6pRcDuv1VXGJnVxlyMJKQN5eXlEHRJ4sknr/GIzZdNYW9vium8bRq7BmDQ1NsKOEuGSU07mij//mdlDh7aocZwE3tYLGHyyTgu2FYo/iZkx0vp9qWqDqp6A8Y78I8Y0bxtV3V5Vv+7IAXYH5PohTCWlealPbNKwJe98DHk3hWDzvX/Hiy++yEHrRy3i7h8n8BZHnVKam/Po27eCGTNmMGHC+BRXFaO5OUxRUYlDAneGkHKVVhJVKG5RGgchY3TgNjHbhJ0SLeb+yRHDW4XiVbOgtJRYLEa4oaHFTtwm8WaMvXg9EGmCoZM/4emnn+bs229nAUYKt1+myaRwJ9KVvFPdgl7zAvAJHAAR2TaT8mnbgQOoah3GC9OHB9p6QDNRn7glcOdCA0UBKCiAL5vhzAMP5JarL+GojfMTVy+QcoyFpZG8f/hhCePHj2fmzJn85jdHeYwyTnerV9dSXl5uEXgfj7JxxBACgQDRaJRAXh4tFtcWM9uqkICl11Y1Khc3cbdJ5CngpQN3S+D2ox8qKaG2trbFhtsJex2hekyUwuFBeOpPF/PUdzN4ap99kFdeIR9D9MkWefB14VmgF7rSi8hDXsmYub4t020nIwJfW9ERUk6u1CdOm2fbdLCw0CxasLT/eEaOHEnZz28gO/WNk3ffcqzlHIBiRIqZMWMuEyZMYNq0LygtLSPuAtMa1dWGwBcuXEjfvgPxlmUNbB14c3OzQ4Xirdg2nvFGheIWrNMRrmyJ2V3Ui8BjeF9dUWkpdXV1nlO4YO55nbXPi8D2NSu59Pzzuenuu5m4wQaMa2ggTOv1NaH1/9Ddrk/qKZAjCxMRKQQmYdwkgph4Sm8kKXsecJh1+pSq3uTIGwU8bLURBY5R1Z8c+YOAR4GFqnq0q93XMB/LNmpVdT9X9zHgnx7DOsIjLSl8As8S2ZJ+puoTpxrFdnYPhUALYMKvD+TVV1/lwAlNUDYEyvoYAqeMeIAqs6TxjBk/cdxx+xOJhB29OR+aON1VV9fSt29fZs6cSVFRCKMh9kbCJKYk0zzHezQRECTlV3E0Ftec2DwfdZdxjgGjA49FowlxUNxhZ2NAYWkptbW1iEebYK5UiL80K6Kw7LFH+fBXvyIwejQyfXqryeUgiUG27Dvr/F8n89LsaPSIF0duPTGvAkRVt7WcDz8TkfXc83QishfG4m1TK+lLEZmuqi9b508CD6jqwyJyDMbseWur7hjgfmB5kjEsdpO6B85W1Vp3ohXtNW1kLXyKiE/+DmSi33RLavY+mfok6NjnB40EHgzBipfu4tHrzmaH9fsY9UlZGZTaC6rZyxubyCl9+5bz179eQ1lZEfX1SQJgO0by6qvPMn36ZKJRaG1RTYt0HQgI777xHz5+54WU+upYDEJ58K9/TGLqe68R1cQ8e28fi8LNV13F3K++SkpAdnoAOOXYY1mzYkXScrYNzc9z5nDxeecRJFFP7nJTohFjlbIoAuvH4OmjjmKb6dMJkfg/8fqfue5UAnqNTjtT5EAHbpk0Hw/8HVqWgJwGHO5R/CTgCVVtVNVG4HHgZKudTTDE/rhV9nFgI8t/Bcxc9/7AzGwv1yZvEdlcRL4QkVetKKu7Z9JORr8nEXnGehsBHCEip2ZSvyciVw+c14RWMnJ35s0pLmYFcQm8wIrXMbgvnL1+Pc/8egF5ZX1YES3kkEu+ZdqMWhKlb0PgJ598PLfddg0jRozgyivv4IILrvYYmTneYot1mTTpaAoL+zNx4lN8//0CvMibQIAhA4vZaOhSDjjkOM78uyWQOBhZHf4XY8ogKnlsuvMeTF1pikU9nlGNGX30ZjvtREF5ORCPMOjcgyHe4JgxbLD55hRZMbvtdDciwKANN2TQ8OGUONJtIrfNCu14iz9tvz19LrqIuRjRq4hERx4vAnffUR9pIjdmhGMwMZKdjoXT8dYpb5Wi3FbAT6omwLC1n2Pnq+pyi/SToUREHheRD0XkOeuFkAwnAgcC71v+LdunKNsKmf7WvlbVhwGsff8M6/cKZKv/dlPpEhH+/tprDD3sMBqw9OCO+CFFRRAqAIqLuPo/q7nt3je59NJviEQKMcRdgNtW4quvvuPKK68kGMx39Wb3bvwcv/76RyZOnMioUWMYPdoZIqq1Id+q+iKOO+44qvo4aNMlLcUstci49Tdh2223pcTqzn4mnSrQmMXQ++yzD1WDB7cpgfcpKuLII4+kmbi0DYkLPthXOXLYMPbbb7+Ess7yUUfa8KFD2WeffVoUSPmOzXbnT0XemfwOejXpZ+ZKXykiUxybM1id/UNd7UirxjvY3cAU5dx5qdrxwg/Alaq6E+Zr4CPLZ8YLc1R1HvGfofdnZBJk+rsJuc4LPUv1IHTmg9OW/tveBzBS4I7HHMOQIUOY/t//0h/LmSRkLFCKi40kHijrw+RFAcZsfSJffvklRx01gWCwD0YCt11PAogEWLOmlqqqIbzxxhvsvff/eYzIshwRYerU2WyxxRasWvULhYW2y71Lx22RdLP0Z9GiRQwtt9Jitl1vNGG1nfpmKC3rS3V1NcUB8zw6pe6EL2WNu93bpOreY4pRv2oV5eXlLVYkXrBfLauWLmXgwIGtVlFzk34EaKiro7i4mHrX/8aek0hF4u0l9F6FzOzAl9sB6qztAa8WXefJJmdS2ah65aW1MqOqXqyqc6zjFzHS/ZFJim8kIocCg0VkP2CddPqwkelvqllEXhKR20XkRUwoiV6HXKtVvParKyu57sYbOffUU1m/sdFMrNnR/QriEng0FOLWD/px/PEn8MQTN3HIITsRtxa3tbXGP/H1199in3324Z133mb77bd29GZr2+Pn06cvYMyYMRQX2/TlDVUlUFjF/PnzGVlVEM9o8dSMq0pW1kN5eTmrV6+mSLyfz2gCoceQQCCl+iQGNKxaRTAYJL9Pn4TFj53SNVZa9dKlVFVVtZC9s4xz8tMm8JKSEhqJT1K6PT2d58nUYqRITwe9gvhzYwduRyUtd6SV4x2tdKlHuWVJ8lK1kw7mAyOT5F0C/Ab4NWadhcszaTij34aqXgPcCfwC3Kmqf8mkfk9Dex+cZFJYW5J4DXDhTTfx5ptvsuqNN7AVIvbiCLYevKgI/j5FOfWCm7jttlu59NJtELHj59nLPsR7fPfdj9h2222JxcIEg16BrOICRlNTAXPmzGH99ZOEubEepmhzhNK+g5k3bx4jB+SbxRwiESOFRyItX77RCNQ0mVCu1dXVFAXiZN1sSeJ2OTeBe0nfTiKPNDTQ1NREUT/vNbad1ijVS5ZQVVXVstJ9y+W42o0CTRaB2+sQRYhbp7j14KngT2i2ASVXBP4DsBITGM/G+sBkj7KTU5SbDIwWkRCAtR+bpJ0EiEiViBztSh4ILExSpVhVD1PVDVX1cDJcZzib39BizMzu946Qrr0a2dzEVCQ+dJddOOCAA7j+nHMYRlyGDgZNPBFbB74mCtNC+zF8+HDqV77G+usPhxa6Dzpq5hGNgmqAKVOmsMsuznkSd+DVAI2NzVRUDGXq1KlssYXzi86iNsck5eKlDQwdNpx58+Yxon+w1QPnnKRc3WQk8Lqa1QQ0rjax986yihU3xSJwL+kbKy2EWWWnuF+/VtK329xwzbJlhEIhivr29fy/OEm80SLwJuIBdSExjG1bVkRe8InbDc1JNEIrbtODwLHQsujMpsDjIrKeiLwtIvan5iRgoogUWrbjh1lpdvC8r4CJVtmJwHdWEL62UAycJyIl1hi2wkyKPp2k/O0iUmCVHYURkNNGplYoFwB3YPQ54zALEPtwIJsJTOfDviIU4s5Jk7j8T39i0OLFLZ/p+VgLHRTEVSlX/6+cK6+5gT9feQGXXLwVcTsJW31iWlYVPv30c3bbbTdefvll9tzz1ySqTRLx1Vez2WyzzZg2bRobbeTx5efQc8/7ZQ0jR45k8cKfGVBi6UUizRBpJmpJ1LYUXt1oJPDGNdUtErdzc+rBbQncSwfuJPQocQIv7dfPU3GZYAve1MTq1aspq6pqcaN3bk401tcnSOB2BERb7+22QvGRJXIngQNcBYiIfIax5Z6oqouBvhjnnnwAVX0NM8H4sbU95LABB6POOEpEPgSOAQ62M0QkT0TeA44G9hKR9xzrFSwGngPeFJEPgLuAA1OEvH4HuFlETsGEzP46nYu0kaktd6mq7iYiF6nqu5n67ftIjTBwxAUXUFNTw9uTJrEuiUSRZ8XQLgjBt6tgzAFXMnXqVPbcdBF9+q9H3NXHjp4SJ+lXXnmNiy++lNdff5Xy8nISXU+ijmOYOnUmv/71Cbz99kvWBKZts0F8b5H4/IW1rPurkTSsXoSoQiTaog+xnzl7stKWwMO1q2m2NC1OqdsuF23porUOHBKJHOtKbQJvcOQ5VSd22RCwdOlSyquqaJ7d+plyknhTXR2FhYVEAgEisVjLHXVOXkJrKdyJgKPvrkJ3GENqqHnT56IlY953tEf6Z8BQV9otwC1J2pkL/F+SvCiwa4r+Lyd9XfZkzIvlTOA6MlycJlPBIXHxxOQxvX24kEplAuaGhsaO5aKLL+bck09meCzW8okexJB3vrUwguTB35ZtxPHHn8CTky7jkN+OhTx3sFMju6uaVpYsWcbKlSvZcMP18FKbOEfz/ffzGTt2LAUF9Xhqih0S0bwFdYwYMYJY/cK4/tsicVuytvXcqy0deNiSwGOxRAk9QQ+Otw7cTeRKXALvU1GRVHVi720Crxg4sOW+ewW+imEIHCDfmsi0Jz6dkrfXyj++RJ4hciuB9zS8DszD6OBrMYvlpI1MJfCo5edfLCJbA19kWL/Ho6MezChQXlbGFVdcQcO0afQhbuUgkBA/e0UzSEkFl116CZceGEAKbPJOVJ0Y88E8fvzxR8aPH8+LL77I/vu7QzK0RkNDPj/88AMbbGALLM4pPjCLORi2XbgsSmVlJYXNiyESbEnXSCxBPRKNwJqwkcAja1YTjRpSd5oStujBrV7twFfJpG9bMg8CK1eupKxfPxY5rsPtRo91d5YuXUpFVRWLXdet1r2OWuXCDQ3EYjHyi4sJr1nTEo/dSd55rro+skGvXtDhalW1VdH/FRHvyZkkyDQa4ZUisgdmibOvVPXNTOp3N3QXKcn+JG+cNo3Xpk1jIHEyaDHuC8S3oaUwaeT7/C/yLeuP3djoVVqFV4qrT8rKyigsLGD69G8544yTqK1dRUlJoRWW26lCMaQ5alQp1113ESecsAM//7yU4cOdq/LEbbyJxVhvRCFnH/0r+uTVs7K6gAqM/tvLB6MyBGcesh818+ewqgQKXJOXLXpwayxXnH46q375pcU1zWsSM4qJX3LbNdegdXW4vSWcqhbbkuSaiy+mqbqaSlr/Bmwijll5++y0E5FVqxL6dlvEO4ncLYlnQkvdX9XRgeilBO4gbxup7XZdSJvDrIWF9wBGYBTtb2fS0dqKtiwNkjlyuPNtSbuvR5k84uRtLwZcHIJfb1xiFpoMGPKORJR33vkC90d8ZWUVZ511Nvfddx8AZ555ISeddB6ffuptFXXZZUfy0ENnMmvWQu67by7HHvtPWmRftVQkkQhEmjnpN1Xcekg9c3Qb/vTtPrwzvSFBfRKJxM0Ed+0P4Z++5bxHXuSZVaG4zjuSqA+PYhZXOPnSS6kYOjRhwtItfa8GNj7tNDbcdluGLlzYctXOoFY25peVcfHjj1P9ww9ULV3aUjbmUfansjJu/Pe/+fHzzxkdDic4Dzlfrs7/qS+BZwnNjRVKT4KIPGHtfxKRH63tJ+DeTNpJi8AtsX4y8GdMsJU/A1NEpDyjUacJEdldRO4VkatE5EqP/EIRuVtELhGRh6yoYz0Gya2v4+dOiW4B8E5EEqK05uVZK9ZYzB6LwZlnPsyUKWu4995HHbUTe/vmm2/ZdNMt2HnnXQmHW/suqiMg1fvv/8Bll11GebntgOtYjcfaRJWvf1jD3vv/joEDB7LZIGmx/3ZalNjneQNHs9lmm1EeDSfYgLdYnhCXgjfZZBPKyss9dd82ioB+/fqxyy67sMaVB3HSjQJaU8Mee+zB8I03bslza/ix+tGaGn71q18xcL31EkwT7e8ab/sdH1mj9+nAr7L2t6rqGGsbDVyWSSPpSuAXA39S1e1U9VBV3c5KuySTztKBiBRj7DHPUdWrgI1F5FeuYmcD81X1OuA2rOhjuURHqVfcE5fuyTN3/pcFBWx1yy0c9rdH+NpadLdl4Rtr6TIV4YILXuXIIy9m5MiRFgl7X8E999zPcccdx+uvv8rOOzvtwQ2d2S+JJUuWM3LkZrzwwgscdNDGLflxk5EoxBSJRnnvy9Xstttu/PLdh/QLRVok6mYPNUrJsLH8+OOPVAZbmw7GYnF7FwUikQgSDCaVvm39908zZzJhwgTs2Jxu6dsm9TLgo48+YtMdd2yxwUmGMmDKlCmsu2Xy2Ppek5fQ+s4n+wrzYUEtK5R0trUEVqREVPUuEeknIluKSD9VvS+TdtL9PYllN+kcwBvQ5nOQDbYD5qmqHUPoY2BfV5l9gU+tcXwDbCIiZR0wlg5Hss/uGmDxppvy0NSpbLnVVjz85ysZ6o5EY+GK66ey774XsXz5cubO/ZTTTjsaMNK0U6KeNWsWY8aM5bXXXuPww//Qas3JOGI888x7TJw4kbfeeontthtLCy062FabmyGmfLe8H1VVVVTUfAExh/t8JHGisjkKFcPH8MMPP1AZ8J68dEq6kUiEvKCZpnGqUXAcK/DzjBmsu+66rMF74tJGGfDxRx+xw447tpB96ys36ANMnjyZzbfaijqPtpwPji+N5wC9TwIHQESOAL4DHgKmi8iRmdRPl8CTKdbrMuksTVRhVKA2amgdBSydMojIiXbUslRRa7oKXvSpwIK8PHa69FJe+fBD/vnII1yz667svfgn+geNsOLk3Rsem882O18MwP/+9yQXXXR0Yh+OwnfddTcnn3wyzz77X/bYw/6o8VYkfPHFEiorKxk0qMGSyu0JzKgxF7TUJ9GmZgoGb82nn37KTiMiCe7zTi/LSAQW1cPosUYCrxBD2k4rFVv6tkk6EokQCAZbBa9yj3bJ7NlUVFQQqqhwXIFB1FEvD/jio4/YaaedWE3qSEYh4KvJk9lqq61ahaXzkWMoaEzT2tZC/AYYpaobA6OBAzKpnK4Vyp4i4mXzvS1wbSYdpoGlGAHIRhmtg8ikUwYrUtkDAAGRDv/vp/s29FKdgLGmiK2zDv949FFKSkrYa/vtGfLNN2yG0ZbYwrQtjNz3Vi2j9rueyspK/v3vK7nppjOI01csgbznzZvHgAFVfPTRR/z+97/B+3aYtB9/XMBGG23Hv/71Lw45ZCsSSN7hRk8syldz1rD1djvz3nvvcu7oPGiOJdh0OycxFzfC5mPG8OKz/2FToNHhuOOUvm3JOhqNkhcMppS+AYKNjcyfP58h664Ln37aUsYLc7/4gr59+9JvzBjkxx+TlDKYOWUKm2yyCXWhEAGP+QIbznF59Z0LuXFtt1BZC4XrdDHNEXe8UUQmA4jIUFVd0FbldDknjJG23VtzVkNOjU+BkXZ8AGAH4GURqXCoSV7GqFoQkY0wJo01HTCWDoNT+lZgFbDLySfz/tSpfPD++0zcaisGf/MN/bDoM2YI3Cbxh78Uina9ngkTJvDPv53FjTfuaUnm8ek+4zBmVCl33nknp59+Ok888Ti/+c0+HiOKT/c9+eRbHHrooXz22dtstNFw4gTu/JyNotEo731ZzW677cbC7z+iPD+aIHG7JzEXNsDYsWNZNffH1gGsiBOzUwcecOjA3QRv1ykGZs6cycgJE5JK3zaKmpv5/PPP2WTHHVNK4AC1P//MqlWrGLLRRknvVqa803t5KjnU+rBLZ1sLsaGIXC0iR4vI1RjuOxK4PZ3K6UrgF6pqK5szxxJDOYOq1ltxAe4UkWWYRSTeFpEbMZHGrsfEY7lZRC7DxM89Ltfj6Ei4yXsNcPGdd7Lvvvty4D77sPrjjxlJotTllL5fmQd5h97GnjvuyE2X/I5Jd+1IIBAjHrHDqRaJsXjxEgoLC/nqq6/Yd9+9CAQCHuVM2VgsxuzZa4jFYowfX0gC/dmrLzgsUL5f3o/+/fvTv+4riEUTiDthEjMKS5oDjBo1ijXzfmghdVv6jmCkAWdMElsH7p64xDXqEmDGjBmss+66zMYKdpEE9kTm9jvtxNePPpqwMg8kSjTOicyVU6e2/D/cEncy+GSdHhTz++ilGIxZmm0k5qe1CKNKqUincloE7kXeVno60bkyhuUg9KYr7ULHcQNwWkf03dFQ17FgdEH33347t156KQNqaxN0Q2BRrUWK4Qi8HJ3A3XvsyV9P24e/nT8YQ3G2BjkCxCzpOwrkc+edd3L22WdzwQUX8PDD9+P9kW+eoK+/nsNOO+3Ok08+ycSJ2ztHgB3/xA5kFQlHKBy0FZ988gk7j4q1TGC69d8ttt4VQxARogsXEMlzTHRao1biJG5L4HnBYEvo1/jVOb8XTOzFWTNn8us99uBfQD+8pe8Yhtwnf/ghN955JzdAKwJ3oi9mInOLrbbiufvvxw5Wa4ebdf4fs5XIfdAytdJLcaZliJEAEdkgncq+VVMnQh17J7EIoD/+SEVtbYJFg1NdYC9aojG4sHIGNx+2EffvvYqQxCBs014TNp2IKCJCLBbh22+/5YUXXmDXXXciLy/P0XprPPnk66y77rpMn/4Zo0Z5rCDlUFZOm1XDhA024b333mXHMXktBO4kceeEZvHwscydO5dyjcXNDIlL3U4SFywzwry8lteTk7ztOva9nD9zJuuuuy71tLYFd+P7zz5j+PDh4Aop63SqysOQ/TeTJ7PxxhuzkviXk/1F4Jxcdd9Rn8gzQy81QsGLvK3079KpL5piJfG1CQERdX9ap3p7teVh6XWeyibYHfzI9uILOjY7lqAd0bsYIwWWAQOC0L8Chg2DIUNg1CjoP6YvjB8H49eFMWMgbzwwBBhm1SxFNUQ0qjz77LP89rf7EwrZ0T5sx/ImDD02AE3U1q7ghRfepbw8yD77bExzczWqdYRCTRCtgZoaqK6GmhqWz13Kc2/OZ94vdVzzK5gxexVDQ0rtGlOsttZsa2qhoR5eXCIsDfZFV1WzvkKfZtOzPYIwcQl8LrCkspLmujo2aWgwErl1L9368BVAXSgEeXmUNTRQhndwKltN0wDU5+fTp7nZ0/NVrLIhYFkwSAAYG4kwlPgagtUY06fVmMmgemgJeGWv4GO/cJwvHreknsmEZ3v4q7O4T030tIyw2UDR9/6YXtny25iqqsmN83sZMpbARWRjEfmViAzzF3RIH07p23nTYh6bO5JeDMuPwZJow2FrfddwEzQ2QWODOW6JXB1G1VChiBIMBvjDHw4iFLJsqlU5//zL+fe/n8OpAFBVSkuLmTjxV+y9t3lG7rjjLU488QPeeOP7VtdUWRbk+F/35+rflDJveYS/rDiUE6ZtwpK6ROnb1nfvUqRs2VDN9pdexeKDDmU13uQdxhDj46+9xtCxY1tI3iZF+/UTA37s04db33kHFWFYQwN98CbvxcB2555L8+jRlAPDUpD3gsJCdrv8chYEg4yKRFgvEmEQhtBtPbi9uVU0mRCyjzh6mwQuIgNy0U5GBO4v6JAdvB7qVL9FIXHSzv5UjzTTEkc7HMb8CTcZFUpjE9DAqlXLefLJF4GoQw+e2Nvf//4Iu+66J59/PpNp075tGZXzfSwiLFq0ggULBnDBBVcxZYo7dp+B/QV3z6cB7rjzLoKFpZQHE1Uo4XB8m9UIBx10EFpTgx1pPEycwMPEFUEVFRUUFBcnkLZzsjMC5K1Zw6hRoxi33XatVpV33uciYONNNuGI009v5ZjT6iFoauLQiRPZ6g9/aNGb2+HCIK6+cUrRaxG3dDp6aTTZB0XkRcv6pCzbRjKVwEtVdTdguqq+i/mS9JEEtn7WDlSVDryk8RhmanqexdU2GTY3xlok8EhtHffc8w6XXvouP/+cx5QpUxFxrjVjWps2bRrz5i1h3LhxxGK1bLbZholjdqjUrrnmJS677Aquv/5yzj13R+9rjMZYuSZCdIPDmTlzJts3f4xqog24k8SX9algwoQJ/PzxxzSC59YiiYfD5IVC2N8W9hYlrhPvA7z11lvssPvu1JPcOacM+Nsdd3DcccexprTUM463rfsepsrtt9zCORdcwAKrjL26qC112y8Lt3rE6+XhIzWcv5e2trUFqnoARhDOB54SkX+LyMFiFrVNG5kSuL+gQxZwTn6lglOia8KsHD29vJw1v/sdm917Lz9usFlLnJFwGJrCQGMDb3y0mCPPm8uWW17IBRdcwsyZ37F8+VLilimGvKurq7n55ru58MILufLKy7j22tahbEQEEeGTT75no4324t133+UPfxhIYaG3cZ6IcP8HjZx+5rk8dPdNHLhOXFqKRKA5HB9rOAyDd9iRb7/9lqLVqxNIu4lEEo9hCDwYCrVI5hHiJG5L4yHg/bffZvfdd6eaxJeek0QFWPrFF3z99dfsffTRrf4XzhdsEPjgsccYPHgwQ3/1q4S5CoiTttdEpo8soPRGCRxVXaWqD6rqPsDpwCDgBRE5PN02MiVwe0GH/UTk33SMK32PRSpSSFXenuiqARaHQvTZZRf2/MtfmPT550xeupRzL72Umpoalq9aTVMTNDaa7dtFcMx/Slg+/BLue/DfvPDCczz66Dncccfv2HvvbXEa6KlGOeeci7n11tu48MILuP76P1FQkO8xauMBec89H3LEEUfw4ot/Y//9LUcWZ1DygLm6huYYc/vsQXNzM+OWvUDAIU05VSdNjbCkATbfaSc++OAD+mF+PI2YCcAGEgk8QpzAnZODYVpL45PfeYctttiCaN++nqaG9n4AcM8dd3D6GWew2KkusvbOpTCGNTVx1513crolhdv6b0j8AnCqupJJ3e70tYyHcoJe7MgDgKouVtU7VfXXwFPp1suIwFX1SuBW4AXgflXNtRv9WodkJO78DI8NGsR+55zDva+8woyVK5n0j39QVVXFnbfcwlaDB3P2llsy5+KL2ePHH4lGYWUD3DGrmLfW/TN3Pfc1dXV1XHjizpx46DKuumo3SkrAUKKtMY5x6613cuyxx/Pcc8+x3367MWrUsKRjfvjh1zj55PO59dZbueSS/0PEsQJkwHFFgQCPvreKk844n3vuuJljNjTUZDvx2JtN4rMVdtppJ/734YdmAQviJO7cbGK2Cdwtdbu3hmXL+Pbbb9ls110TLFRw7QPA5Oeeo6CggA333psY3uQtGAugl+67j+23357mjTduWRXJlvCdiik3eXu9PNJBbyV27aUSeDKoatrKokwnMX+jqm+o6s09fTWeXCOZ7jXmUcbe27rxqoED2WHHHXnphRfYcdNN2WX0aO4/8USWPvMMo1esYBTGfK0J+LBW+GSbw/nrW7PYcaedOeGAHVnvh8u4/+R+jKwQoBaRWgyBNwDNfPDBe6iGKCkpYdGiuey7756JY3fovaur1/C//61kxIgRNDZ+y4QJ9rJqARCH9B0IEFXl0xUbMmzYMEpmP0nfwrjOO2qpT5oaob4e6sKwpKSEzTffnBkffkg9hrzdkrdbAs93EHiyrQh4++232cXSg3sRqo2qaJR777mHU886ixVWmpu87YXpRlVX87e//Y1jLriAX6zgW83ENy9VTSpJu5fwT1Zo8Q9rY/ORiExVKJeJyE0iskmHjGYtg5OwnV9/NnHbUvjPX33FWQcdxLuTJqFz5jAIQ0q2JOk0n6sPhDj48CO59Jwz+Oj83bh3w6/YcWgE6uugtg5ttGnR2HUvXryAhx/+L8cddxy33nojl19+IW4qsQlcVfnrX//JlVdezdVXX8HFFx9glXAs0+ZQo7zw8XIOPeYs7r3nLk7bVhPiXDnVJ+GwIeXB227LTz/9RHjJEk/SdpKyEpfAI7QmbacEHgDef+stfvWrX7XowZOpNELAC3/7GzvssAOl66/firzt8wJMeMv/3H47Bx10ED8OG94yaemU8lPZeLcFn9ANeqkVCmAipranfqYEfhRwBbCNiNwjIr9vT+drE9yqEi/ChjhpOyfXbKedZOaDNlGEgc0aw9x74B7sOf1Zfj3AmIBTWws1Nejq1ciaNZjFrWuBei644Fr++tdrOeecs7jllqvIy7NHEg90ZZsPfvvtHPr334AZM2aw885VlJWVkrjGppjl24JBCObz7Bdl7LTTTiyf/AhD+5qfkm3zbevAGxuhPmycbTbbeWc+/PBDSog7vjjJ20nitgQesiRwt9rEPak59cMPGTt2LH2GDvUkVRsBoHTVKv75z39yzJlnUksieTsdqwqBMfPn89///pf9zj6bJSTqv5NZoODRr48UUGMim862FuJwa3Ux95oHaSEbKxTbZ3sr4ORsOl1boI7NDbeU7VSbOCXzZHpT28rB1mK3EEZE2TUf8qJGNVFfD821YaitayFyQ97GP/DGG0/lzDNP4LjjDmDgwMpWPdlWJ6pmmzPnGx588HaOOMKOF+6UT4OW9J3H8uowdeEgl1xyCafsDPZSPk4LlEjEmnAF5gEDBgzgw3ffRfCWup0EHQMaGxvJCwaT6r6dW7SujhdffJFBY8d6hPNKRB/g4TvuYMXy5awmUfIOQsvq84XAuDx45KYbWbJoEQsxP3ynTXpb+m+fyNuG0qtVKIeo6umAisgDInJGktDdnsjIlV5EvsSY1L4K3KOq0zMdbVchV670bmnavRqLUxK323CvWB5wbUHH3ulWX+jYSjDEUwGUF0JVFQyoguHDYMQI6DeiD6wzFkaMNPs+wzHrT1ei2g+RPhhZ3x6FU86Pu9THYnWsWLGcAQOKiNOsrbGuA62F6mp09Woii5fy2We/sNOAeli+nKbF1SxeTKttObAQQ+I11nXYpO1UQ7hHVWfdzxJHnvveOm25GzB2rSWue2nfW6y9WGULMKoqu0yelVZitdMXKC2GX8IgEfP/rrauoda6K/V4v4SSWcO0R0fenpdBZ71IsnGl37BC9F97pFd2/afXLld6EdkRmIMJzncYhltrgZ9V9Z626mcqgc8ENlXV03oSeWcD9w/eS9K2o4q488W1pduX18PuJDTbpdxWTdgThLW1oDVrWmKUGCnc3uoQsV1hmmlNI4kIBPIZMKACZ+SWpqYoL730HRAECUIwD4J5BAtD7LRZf0ulEiQQgEdmBvh0WdyJx+45innz98WlJqG1NG2PsBhDsF7SrrsOGOLNx/veVQ8axKohQwBD1GVW+27yjgEzRo1itgiFIfNhMShg8pJ9KaRSpbj/xz5aQ6E3mxE+BryHCZW9uaqerqoXY6SvNpEWgTtinlwClIvICGtb68wI21KJuFUgXvlupPPweklmbqKwSbwhZumWLfJebQWOMoGmVhkiX11NnMRtedGm0hjRaDO33/4g117rfMm3DtVVW9vEccc9yY8/jufhh6cAeS2ELQUhQ+ahAgjm88R3MOLEv7Po/85nTn38Htmk59Yfu4nZixDb8nR0E7UXedYBR5xzDnc99hhLAoEEnbcz8mAIIwrd9swzbH7JJSyzQt6Gw3EXf2dMlrZI20ea6IVmhCKyn3X4H1WdoKq3qepqKy+fNH1s0pXAP7f2bwP/AB6xtkPSHnEPgpuENclxJi7y0JqkvT6pvTa3FUYzcQKvr4c6o9Uwf1auonnpclb+uAjzwb+auE7caJ2/+uob/vjHc9huu72JREr45htnoKq4cmfFijqOP/6fXH/9Q8RiUUQsp/JQAYQK0PyQdZzPy99HWLr99Wy//fb89OojjC2JX5PbmsZLgnUTopeDjBdZep07/0cBoBx48MorGThwIPufdx71xGOb5JEYBXIH4OJjjuGCiy7ih212obY57iXqDNrr9QLy+h8nw1rERe2GKkSi6W1rEW4QkWeAESLyjGN7GriZNJ150l3QYWvr8ExVfdFOFxGvtbl6NNoiZCX1Wy9GopY5nTek/TC7I925JUt7QjMMNFpqlFW1MH8hLF8KS78vITpkNH1GbELj0o+5557B0LcYo1wI0dAQ4dpr76a4eCR/+9sj3HrrrUA148ePajWmBQtWcN55j3H33Y9yxx23s+WWdfz2tzsC9WggD0L51hbi0zlNfBA6grMOP4Jz/7ADF41dxqIFkB+MTw7akCTXmKlu2L5X7rZt6RoS9eRDGhs5/rDDeOOjj9jnrbdg2rSWIFW2xYk9z7D1D99x/llnMumJJzlk003ZYunShHgsyYg72XE66O2EvjZJ12liFmZpSC+UAzcAB7bVSLpLqtloCbQiIhsDuwGvZNhGt0G6BJvrmLnJ+nWTv03aAeJS+HyMLB2sqmLIZpsxdrPN2HSzzdh7880ZM2YMc+fOZdq0afzjvg9g+XJDskVFvPPOLP7+9+lceeVfWbp0KSeeeBR/+tMJrL/+OGuB47iN1qxZv3D11f/m/vsf4fLL/8RBBw1kl13Ww0z9GTsNKSiEggKm/9LEP77djOtvvIyTD/4/7t5mLtVLYGUIQiGjOy4gHtHPHUDKvlbn3n3cFpzt2S9Y95YHxL76iqv//GceeOIJDtl8c9ZvaGiRvIuBfgEoK4O8IKz49yO8tvMuXPH441yz554MjcVa1CjJdODJxp0pmXthbeY3e03MXoY7rYCAnhCRgek0kukk5gT7QFW/zrCuDw94SZvJ1CjNwPCjjuI/Cxbw7oIFXH7zzay3/vr877PPuOrMEzh68/7857ixDH7/KJ48qgGqq1n502JOP/1f/PjjcO6550Huvfduvv76ZR577Co22GAMifIifPHFTK6//mXuu+9vnHvuWRxzzAR22WUj4rKqpTEOBpm/tInrHs/n5lvv5rTjDubOvZcyoF+AomIoKobiYigsjFvS2KoKp3UIJCf0tuB8GbitTtwWPvkY6fqNW25h4cKFnHTzzaywxmNPltrjzQvAegH452mnMWjQIHa77DKWEDch9PLAdI57bSbbjkKuzAhFpFBE/iEin4nIFBFJat8iIueJyFRru8CVN0pE3hWRD0XkPREZ7cofJCJviMg/PNqtEJEXROQjaxybu8u0Qd6/wsG1qZCWBC4iZwFnYyYwjya+aEmyT4BehXQkebuMl4rF+cAnU6NEMP+srz78kGMOOIAV33xDVWMjY4OwXiXsUAWDBpltRN868mpX8sy/p/Pmkm245tqn+OSTTzj33EO4+uqJDBs2BnvhMpH4EmvvvTeZZ5+dxm233cEppxzPFVfsx/jxA4krDmwr6SArVoa58C+/8MCD/+bM047jL4fUMZAyWFhH37LVlJfDmhoj0dY3Gk28c7Ua5zU7Az943RMn3OaDbhNMp1rEaZJpW5mso8pFRx3F6199xTuvvELw5ZfpA5QWGQIPBOJWPluGGzjlD3/g5c8/Z/JHH1H3zjsJMcnbcuTxiTx95FCFchXGPHpbERkPfCYi66nqEmchEdkLOAHY1Er6UkSmq6rNaU8CD6jqwyJyDPA0sLVVdwxwP8ZC1gv3AtNU9UoR2Q14XkTWUdWmdC5AVd/GzDe2ibQkcFW9Q1VHAyeq6hhVHa2q41T17HTquyEit4nIX0XkEhF5SkQGOfIuEJFrLO+k3ySpP0pE/m7Vvz8Tw/fugmS/VzcB2JYbNokP//FHyidPZuPGRsYBFRFobowvX2aveMbKlXwzu4HjTzqTP11yDtr0KH//+14MG5YPrCHuC2mm5t5882Nef30Wl19+BSeddBQ33DCR8eNH4jQnBEE1j/p65fTTP+fue57g8ssv5qxDGhmzTiWUlkJZGX3KoG8Z9CmD8nIoK7UXeDNbMYkSuS2Vuwk5HXv5IIlSfdDVrpO87fN1fvmFM08+mXsfeohvq6ooKzFDDwZpsTppDJs7M3jGDM489VQeeuIJlg0a1CZxt+d/31thq1Daa0YoZpb9eODvpl2dBUwDDvcofhLwhKo2qmoj8DiWY6IVKmRTKw1rv5GIbGGd1wD7Y8yq3WOoAP7gGMO7GNllP3fZXCAjFYqq/st5LiJtKtmToE5V/6Sq12Fu8J+s9rYBdlPVyzES/y0i0tej/iRMNMTrgG+Bi7IcR0pk+6Blo8v1mviyLSpsj0ynNUozcTebNRZ5V1fDypVG9b16eYRLt6rlhev35ZaJSzho+woksgLj1L4KWI3qGlSNeeHGG49l8eKfOf30Y7n33jMYOnQgrenTbPPm1VBePoo777yDA/ddzaZbDm0hb8rLoaKCykoYWAWVlcbpqCoElZitAjNLU4YhdKfDkpPQ82lN8s7zQtexs42CJHkl1hga/vUvXn31VX5/6aXML4CiIkPetnWP/YrLA2Y8/jgvvPACR5x1FquIm0S2ZRHj/l/68IaSMyuUMUB/YIYjbTrg5fizVYpyWwE/qWoYwNrPsfNVdblF+l7YHGhS1flpjKHdyGgS09LNXId5BgTzDD6baaeqepnjNICxcwPzlvrUKhMRke+BXTDha+0x5GMmTydbSR8DfwMu9xjviUC7gsV4Id3Jz7bgdKV3qlWcfbjVKE4it0OyhiJxm/DCQiOFFxbCsNIa/rp7JTStjk9oDgqaqILkAWqpUJSBA/ty9dWH0a9fGaWlhTgnNe0ICqoBVGH99UdwzjkxZs2az66/Hmc6LA5DaR+zQlB9HUW1tVTUh6mtjbvUR5dCLJI4sWi/FrysO9yk5zU56ZbCvSTufOvcJvEiYLMgPHXuaVTkRdi9j/kxN1v3sT6WGKelH3Dj6adDJEIxcQ/SVNYmPmFnAM1IhVIpIlMc5w+o6gPWsT3xt9qRXw2s79HOQI9yVUny3Pmp0J66GSNTK5SJwJ4YUrwVOK89nYtIObAHcJCVVAU4jZJraH3hlUCDxmMAeJUBwPrHPgDGld6rTK7IOB249eBKcgsXLx24bY1iE3gehmDygGCtsfoIBuP7wkKoCixHCosSoghSZTo2S66BTW3Dh/ezenGuM2NTpYmFYt6fyvjxwxg/vgJoMB1FIlDaZK3RadZ+6x9bRCQSIRCABQ0wtR42bIDCprjlh22h7lwb03bsZ/BgmhoaKKquTjANbP1NAMv79CE/FKLfihWUk0jcNpEXW1t5uRnyIQV1BAIQjUJ9g7Gnrw8nLjBh24BXhsMtsVu8VCg+2ocMCHx5Gq70Xk7T6ZRrKy9dg7T21M0ImRL4TFVdJSJBVW0WkX7pVhSRkzB2jbWq+ntLNXIPcKyqrrSKLcUYC9gos9KcWA4UiYhYJO5VptOR6kWQ7iSnE16TmTaJN2H+cbYapcWTsDZO3M59/8IlhrhDBcZzEtAqQcTutdTqwZZjW//WTIgLIzerirXIQ74ZVUGhCRVXWmKFjYuYxTsjEarCi1CFJxduzLWvPcBd+27LVkUQqo5Lw8UkBreKYESYi266iSlTpvDR7be33BN7bxO4LcVvdOih7Lv//tzym99QjLf6pCBodPNlZeYlFwhAUxM0NJitvtEOAZa4Pqf9wnSaDCZTnbSlPvHiqd7+AsihGaHNA+VJjt1lyx3n5cCyJHmp2vFq1632LQe+S6NuxsiUwHcRkalAoYj8DbMyfVpQ1fsxM7eISCVwO3Chqi4QkYNU9T8Yq5YrrDL5wHrAB9b5cGCh9eJ4F6On+h/Gea7bWcO0Rdpe1ijufHtvE7dzYjNMaxILhRMl8YJCQ1LBYJi+gSVx8lY1T82AKOTFiNOmTXXOn0WsJV64UbfYvSmGwK1p1sIiQ9qlESsUYRRiMSQWZVBwOZt+8T1jxowhvMl2DFv6KcXFVuiWeu/AUAoUFRURbmykwnVv7O8B50LDw0eOZN68eQyAFgncfh3lFUBlUdxMMBQy7YTDhsDr66Guzui93SsEJXPeSSV993ZCzgY5ijT4AyamyLrEyXZ9vH1VJlvlcJSb7MgbLSIhVQ2LSAgY68hPhS8w/DhcVX92tP2PTC4kXWRK4Idgfp+fYWZ7b8qy3zesvh+3wqyswcQE+MyyvbwWo3o8T1WrrTpPAxcCH2Fmi6+wbDxHAOdmOY42kWsVi7s9N5E74VS32ERuy8bNJCo3bFoN1BvyNsRtyprzRvoEjCWVgLXKRBSqYpBvU5MdOiofmx5VlebmMFOmzOCtt77gyCN3YdSovo7eLfk/L2oYMhaLzzjZ38XBfE7efiGPPfYYux9yNI3//JRhxWbe0469VV8fl3zDGBVGUVER4YYG+pMYKAziHp55mCnZkSNH8tWXXzKcuL47FIKpoTxGHn8Gbzz3ONvLMtYrMZedEIqgzrxEbLVJHd7knY307aNt2As6tLsd1ZiIPAgcC3woIuMw1iR/FJH1gLuBPVQ1ijGEuFVEbA47DDjfaudLEfkKozJ+xNp/p6pT0xjDChH5lzWGP4vILhhZokOEzIwIXFWdAVbutCYJW5nSpNFOK8N2R57nS0FVt3ccz8XcoC5DMqk5WzVKOlJ4EFgzfDj9+vZl4YIFhFetohSjc+qHQyqvNQSeHzJEvzoGM1cDv9RRH1rAimANK5jPqlghjcF+rLdhX049ZzerpWZUC/n22yW8+eY3fP/9CoqKhrHjjjvxu9/tyR13XEJ9fYySEqFv3yD9+xfRv3++2Sry6F8M/QvyKC0pQRwRiMZtAF89+yh3PPkelzxwNqeMbKBPmZHCyy0rmvr6+ILN1TEoLCwk0tjIQLxfbvaLqxoYNWoUbz7/POtbErYtaa/Iq+CY3XfnL9dey3/+8x8ee/Aeyqd/xqYBaK6Pvzhs1Yn9EnG7zjvVWJlI3z6RpwFrQYcc4Spgkoh8hnlkJqrqYhEZhXGOyQeiqvqaiGyAMYIAeMhhAw5wKPCQiByP+cw82M4Q8yn6NpjVDkXkPeBIh+XJqcA/ROQjzE/0tymsVtqFdB15VmKek4RkjP75gVYVehg6aiIzXdJORwoHQx57//73XHrppVRWVtLQ0MDChQtZsGBBy7ZwwQLmLFjA9x9+yJDxm7PR5pvTv39/+vfvT0VFBcP692cT67y8vJxAIMDZJ+/F3K9m89bny5jyRQMiI9lyy+056KDLzSIMH37Iu+++y+uvP4NIA3/5y2MMGjSI2tpaVqxYwYoVK1i+fDmzZq9oOV+xYimjSmZwzkFDkGAQgnnsNOQXfvjhB3ST3zGkz+OU10LfvmbysK4O1qyJqzR+qTESeLSpgSHFcQktFsNeO6JlTnZ2s5HAaxfMo6rK8v4sNC+w82QZn5+9H/cXjWbHw0/mwf++xPz583ngnnv4/sknWSdWTwgjfbvVJ6kiJ6ajSvH6f6eTlk2ZngzN4WINFlEe7ZH+GTDUlXYLcEuSduYC/5ckLwrsmmIMK4HfpDnkdiGtBR1E5DBVfSLd9O4IrwUdEvKzyPNKd6cF2jh2752rw7gtLsCQTBQIFhRQOWQIlUOGUDV0KFVDhzJk6FCGWtvNl13GyE02YaONNmLFihWsXLGCuuoVhNesQBtWkBdeSVDqCRWECRbl01T5W0444QQmTJjA5MmTefvtt5k372sqKhrZddcJ7Lzz+pSWFvPpp7N4/vkp1NVFaG4OopqPSAF5ecUUFJTSr18F/fv3JxQK0bzsYU753cgWL6M1S6v586c7sdc++7Pmn/vxmxFN1NTEnZBsAm9shNcWwhHPTuOWyy5iz5lvEGk25G0HNrbJG4FnVgd5bnEjB00YwjGlS8kLmkBaQMvqQI2NMHc1vFtfyND9D+aE005j/PjxPPzww/znvvtg9mz6Ya8kGpfAsyHxZKqUriLwznwBZLOgwzrFojeNT6/s775auxZ0aC/SksCdJG0FsRqAUZ082UHj6lZIJklnKrmnUrvYeztaX7K6tvY52NRE6KefqP/pJxZgZmy+I/6ZXwxUf/AB3+dDVR/YoBwq+kFlf6iogMoR0L8SqKxkTVEJh/13Dnfe8mcK5Qd23L6Kk44cR9XITayWQhgqq2O77Yaz3XbDrbU0nYu+RVGtp6amlhWLv2PFomqGbzseiuPicp9gHvVzX2PHHW/gtNtGccB2qynvU0N5XT119dBQH9dN92kwErhEGxg5wpCwrZFxWkSKQHHNMMLhMOXRpVRU0OJN2RSG5nBcLVPUBLs1NbLq0Ue54tFHiW25JRNPPZXXp01jx/HjKVq4kAYS45e3Rdz2AsfJ/t8+2kaudOC9EZk68lwA7IMJivcP4Aw6yAuypyLVJGWqcpok3YZNFE5Vix1DxDGViB08tgDLy7EZpB5qxARpEjGfrGAm8wbEllPaN8wDe4YZPKoS7dcfykoQ6mHFUjPTWFBM3KramBkar2V7JMZuXETo27eQvqX9GTOi1Dj1NFp24pad2B92K+Wll15i2BZ/YKk8xcABBVBSR2l9A6X1dTRZZFu4wOjA82hk0CBobjbj9XrQC4uNBcowyx0+HDaE3xyGuvr4ykX19cZUsBnjslc3ZQqPH3ssN5x5Jn1qaxNMBeOvpfTUJl7SNx75PjzQO6MR5gSZWqGUqupuInKRqr4rItt2yKi6AB2lB0+nn3R14ba8G6Q1gbsRxbF6UKPxMtQYLat7h8PG9rmpCcpraxjUvwGdV4fU1qJ9+kC/NcazsrTU2HcXFsVnBsU24LNfHUni8wUCEMw3rBsqgFiMXbYdwgl3PcGlV97EYzc8ynl794PCAihugnAZBeEmChqbyJvRSFFRESEaGDwEopE4MTtXvW+OQEn/UcybN4+KoCHphsZEab6+3hB5fcSoSOyJynqM3c3w2lqaMHlOiduWwLX1lbXAl77bD3tBBx+ZI1MCt0PX2b/bHhdEqjOQrRSeqr4Nm+ibiQd3chJ5gLgTvN2GArEIaI15WOwtGjVbOAzNzc30bVhBXnMEqW8wLFhWZ8i7tI8xESwuMkRcWGAtq5Yfn1FsGaz9lrB1Hg6xORAgEMpnZNH3PDDpbvpRgBaEkICYthobWozYGyTGRx99RKSxmtJSI03nBc0+HI43uaYRFjZX8/rrr7MRcUuWFonb3mKJxG2bC9rOUDZhu80G7Zdme6TvZPBJPg5fAs8OmRJ4VEReA4pFZGuM0XqvQK704G217yWFe5G5LV07Q7E61+q0t5gzLQqx6rgE29hoyLCuzgrgtAbKa1dTVrYaKetjXBZLrVB9xSWGwAuLDIGHCggTIFSYbxma58UHYpN2i014hE8nL2ab9UoJxGKc+Yeh5MemUZI3CG1yWVdZzDyyHzz/yO1IUz1Ra7y2FG7rtBsbQZth4vLnmTcTSoJQ3QAr6kCs/DWNZkKyjridt9NF3nbht70to47NaTboXuLNPk9nDVSfqFNDM4uF4sOBTO3Ar7ScZzYGvoJWTnI+0oCbpNt6MXj9tm1J2/aFtAnbJhR7vU6bHhMWIqg1D0xT2OiWGxrM3ibFNbVQXr+G0lqHGqW4yJB4YQFaUMgN/13J/34K8c/Lh1LSpyBO4rZEbj2RGo1w7o2zqRj7O778+t+csv8gygsUwjErZoqlBwk3tZA94TCRcJSLr7+f2/50Mg2NHyRMRjq3+jooicDoADTWwwsr8tAd9qDfq68yQOO23TZhO00E7dV1kjnreJG2DeecRbL/US7QW3gtV2aEvQ3p2oGHMHFMlqrqG8AbVijZqzAekmsFspWmU+m1M2nbLYU7JfAESZr4hKZTIrf13nY/Nmm760UwTixNYSPVNjUZArfVDcXFxja7pBT6lK6huHgN+cXBFgL/YXWQH8N7cvv9V3HaKXvw91PKyCsIxV1AAQLmVSIxpShWy8EHH8yD9y/jpxkfMLoyL07WTY1ouBlpbm4JgkVjA4FwhGg0ikqQuloXaTscfurqzFxpXS18tQa2+/MN7PZ//8dpb7zJFpFIgm13y3qieBO3Td5ezjrOtTyd9xISSTZd6buziLknvAB8K5Tska4E/hjGlqzCiklyCsaD+egOGle3RC4nOtORwpOVtdUnzklN594u1zBmDAOGDGH6xx8zUpVyEgkqGoFotZHAHw2XMbC5hq0qYL2KOJFXF0OfUigsjFBcvJrCQlinMMCejY/x6qtbc+Y1T3DhZftw8+8LkPygi8SNWuWKXwc45sITeeCxlzj5D9vxz0OqCVgG2hJpNqRtSePRcMwI5U0QjUaJkUd1tTeBO00Pf6mF+oMmcuRRR/GbLbdkg0iEGuIu8bbaxKnndnpZJjMVtF+MNnG3TAxD0glMXGU7JAzd2gRfhZI10iXwhfbqOyLyMXCRqn7QYaNaS5CNFB7FWoCX5NK4sz3n3knizcCoTTbh3nvvJRwO89RTT/HKk0+y+ssvGU18pZ9moKYJTnjgLvbYYw/eeOMNXnjrDZZ88iZjI0vZrgpG9DMGKC1rXBbG+L/+Me5/5nQWDXue3c94mNvu+g3n7hSIW55AC5EXAmePXsADD9zPSZfey3237c5pWypEIi367XA4bmXS2Ah1awyBq+axbLkxBWxojO9t4m5oMHruKZtsynP3389hBxzA6HnzWiRvZ5haL4k7mdTtVp24pW2bmJ1zb26duE/i6UFz60rfq5CuQLnScfyqTd7WAg9rFdoSBHL5SZzs09tNzO68mCvPK1aHAHOefZbthg3j1OOOo7KykufefZcnv/+ejS6/nBnrrMPXmPV5osADxx7L8QccwA9z5nDMSafw3+9/4ejnpjHv9zdwV+D/uO7LAp7/FmbPg19+gfnz4YABUZ6//BCqqqooPeAeHnw3zLL5jayav4bVv6xh9fzVrJ6/mjW/rGaC1LHk2cupqqpi7uhT+Py7CIsWwkJrW7wYFi2OH9esNgQeJY8li2HpUli2NL5fvtysPrSiAd6q6M/fn32Wq6+6isg77xAmPmnpVKGESU7mycjbfa/d/w9xbJA4F+EOwOUjOXK1qHFvQ7qu9B8TD/qyNSaMK8C2qrpzB40tp2jLlT6hbJb56bjWe6Wl42LvtW9rs80M8zDqgzWhEJvvvTcHTZzI/vvvz/Tp03nqySf54K672KK5mRCGzBcBK/r2Zdj/7cYOv96DPffck0GDBvHWW2/x8PEHcNxIJZRvmYQH4fYlg7j9pc948P5JjP/kejYfYJyGINHKsKEZ/rJqGx5+7k2O/NXmXDZoDrFoa9vucBg+WAa/e/oT7r/5Bnb47Hmawkb6ttertAl5Sl4eV772GkuXLuWuP/6RQda1Ol3g3eFgk0ndTvWTm7iTvUy9jp3SejpWKsmQC61CZ2smsnGlH5EvelF5emVPX+670juRrgrFFmoA3nWk95gPn854eaczmZlpW6n2XuWdcVMijr5DwMBwmMXPP88Nzz/PVaWlbPPb37LTbrvxXnMzq4mvXjMEGLl6NaFnn+PbZ5/jpQBExo6lfKNNGN6sLFmcuPLPIYHFnHLQPjz+5kecf+Y8mr95khElpl+XUQo713zOpPvu47xbH+bBP+7CnuWxlnglzZb3ZCQC1ZYKBcljyVLLUIU4OYeBucBvr7uOyspKzvnNbxiOCQvrpSbxIm0vyZsUx5Ce/tue8BRXWV8STwJfB5410iXwC1W1VTBzxyrNaxUymVTMRfvJJjQzIXG3i30yIgfjat+3tpafH3+c+x5/nEqMRGt7HoYxRN6IIfUNY5A3+wfyZv9AMA9WNFl+PEHIt8zAf63TOfYPv+PRZ1/giN8uYI+5H9DP9ckTjRm702dvvYL99p9Gn4POZPrjt1MZMAQdsdzlIxEjZUejUQjksTpCyzJrtmS9Ehh46KEcc+yx7LHllgxqaEi4hkyJO5mqxC1lO1UjuMong1edzkBP4kSfwLNDusGsPFeiSCfAeXdCe4k32z4ylcxTlW9LEnfCJm57c7cTxEQAj2L0xbZtdD7xtSYbiC8UnI9Zt6G+wbGABJCXZ1Qm63z8HqcdfzzR+gaWL4NGiQeccoaD3VqbOPfoo7ji9ju44M472VONFG7HHolgVviIRqOQl0c1cQK3rUiWAusOHcqRhx5K0dy5njrtTFQl2ahLUgUec0NS5PV2+K702SMtHfjaABFRm3DSRS514V7p6ejCkx0n04k7j9PRjzv15M5V4u2VMfNxkLejjLjSbPvzZUB/DOkHaK02iGFIeCHwswibqrZM/DmDR80H1j3rLKZ+8gnDJk9OCDBlu72vtvYFjnqZSNy0cUwbx17nydLSycumXEe3kSmy0YEPyxM9qzC9shfW+zpwJzJ1pe/x6CopPJ1ybalSnMfJJHG3VO7UjTs3d32bzCMkrjfZjL0efeIyZgEMaRaRGL88hJGenUFzFJgNRAMBJliieBEwVpUa4p6kNrlHMWqeWXfcQRlxvXYMI31D3GRSiLvAt6UeSVfi7mry7nXwdeBZo9cReC6R6cugPeUzIfFM2ofWpG6Tcx5xNYxN5DZZLxg4kL5jx1LyySdY85Ut+faxrfv9ATjh4YeZ8vnnLJ00iRBxG2qnu7r7XInbq9vHzuiAzY7jdIkbj2Ovfapjr/Puhu4+Pjd8E8Hs0NHCaEqIyJ9EZLkr7XARuUVEbhSRk5LUqxCRB0TkYhH5u4gM7IjxtechyPZTOdW513GqvXvzmtxz20Y7N3f8EHurBkZtsw3/eflllm6wAdXE15WstbY1QI113B948d//5pLLLmN6YSG1jrwajCrEttuud7TlDDxlB6JyusM7nXRSWZ6kInKvfapjr/O20jNFTyPf9kKteODpbD4S0WUELiK74gqGJSLDMCtDn6+qFwLHWytLu3Et8JaqXg88B9ycSd+5fECylXpTpeWCxL0kTi8y9yJ22xIlGaED/PDCC1x37bU888or/Dx4MLXEydfebEJuBma+/DLz589nz1NOaSH8BuLE7IzT7YwW2Ehr88FUhO01YZmuVN7W/U6Ftsr1NlLOFHbwyrY2H4noEgK3JOZDgLtcWXsCUzU+s/opsLdHE/taeWAcjPZN0s+JIjJFRKZkO9aOkMI7msRTEbeb1NzHbseXZNJ5IfD8TTfxyiuv8MiLLzK/pKSVpG6Tch3QD7j6ssu4+JJLWFhamkDUdkzuMK2lfXt9SvdYUhF3W5OW6dxP93GqNB/tgwKRWHqbj0R0OoGLWYvrWuBSj+wqzNe3jRorLVW5GqCfiLTS56vqA6q6pdesdWdJ4bl8AaRD4umoUtxEl8xu2i2Ru0m0D3D96aezdOlSbnniCX4JBBLI114U2C4/6513+Pbbbzn4zDNZQyJp22oSN0knW1jYmZbOBKb7nnjdr2zJO5fSdy5+lz2R59r6SnL/f3wYdBqBi8hJ1mIQb2Cev5MwUQ2LLF32OIx5bx9HtTIrzQ1nuTJglapGPMrlBB3xw0mXIDIlcedxpoTulsSTbU5yLYtGOeWQQxg5ciTn3HYby2hN9nbQrFLgmssu47zzz6emvDyBqJtdbdvHbTnmREmPuJORufvY6zxZWqp0H+nDHXfdJ/D00WkErqr3q+peqrq7qp5s6a/vAxpU9XpVnQ28DmwhIrYt6XbAqwAiMlhEbGvRl608gB2s84zR1VJ4R5B4KukyVZpb8k4loTuJNgoUr1nDwfvuy0EHHcRvLOnarXKxbbxnffIJn376Kceedx510CaJJyNvdezbIoC27lmye5sqLV34pJMefALPDl3myCMi6wAnY6Tw64DbVLVORA4HtsR63lX1fqv808BzqvqkiFQANwDzgLHAxaq6pI3+1MtmsjMce9rK90pPJy0bR6B0y6Vb34YC5Ztswpvvv89RRxzBty++iFfwsAhQufnmvPHee2w8Zgyh5csTSNYdayTZl0WqfTbHXufJ0tLJy6RMNmU7o51skI0jzyARPSrNsjfiO/I40es8Md3IJYG3Vaa9npteae05b29b7rQIsM5uu7Fo7lwafvopIcyqE2uAe/79b95++21evO++BGcE2xbcyyGJDNJSHadzniwtnbxMyrSnfEe10R5kS+BHpFn2Zp/AE9DrCRzWPhJPp0x7ydoLAQyJ2671zoUMnMSiQF1hIYHGRgqsNLtcNEkd97HdtheJp6qbznmytEzy0y2TTdnOaCdbZEPgA0X0sDTL3u4TeAJ8T8wMEaN9EwfJ6meS7k6zH9pkbvleZZwPerpenG1du91GOv4WBY2NLXWcZOxcLs4+d4/Xi8DdY3Afe51nkpZJfrplfMTh36/s4BM47SflTNvrCBL3SktG7MnSnOleD1Qm7vrOGNheS4q5V66xJ6nc4xLXubPtjpayvdCdybunkqD9//eROXwCt5AJiadTtrNIHHJD0G0Rd7K8VEi1tJiTjO1ybom6LQk7nfNM0lKlt5XXHvjk5d+DbOETeAciW8m+o8g9nXQ88tvzcLmtSrzCy7rz3Ppy5z7VmNIh82zS28rLppyPROTqvlmmxpOACRh+u1RV30hS9jzAVr8/pao3OfJGAQ9bbUSBY1T1pzTrfokJGWRjjqoe364LSwKfwB3ItRTenjZSkTUeeekQcypSbovQk5VLB16E7IRzdfdkUnqqBRFyodvuKvJeWyYv2wPbpj9HuApjnLGtiIwHPhOR9dxmxiKyF3ACsKmV9KWITFdV26fkSeABVX1YRI4BnsasB5xO3S9V9ejcXVJydGk0wp6OXOhDc00qsXbmpSrjVS6dzQk7LKy7jJu4vVQpXjrw7kLePrJHrjwxrTAdxwN/B1DVWcA04HCP4icBT6hqo6o2Ao9j/FIQkU0w5Py4VfZxYCPHEpJJ63Y2fAJ3oSMe2PaSeEeSdVv9ZvMwparjpRJx5nmpQTIh6VySdybwpe/2IRcEDozBRC+e4UibjnEMdGOrFOW2An5S1TCAtZ/jyk/Vx2AReVZEPhKRJyx1TIfAJ/B2Ilef1x1N8ulK3pmQQZYPGTESpfB0x+x+CbjzkvWVahypkKv/ba7qrM3IgMAr7Qij1naioxl7XYDVjrRqvAPiDUxRzp3XVr67jznAyaq6I/AJ8KGIFHuMod3wdeAeiNExK+20Va4z8mmjjLOcFzKZI0iFdHXubnvwTPrpruTtIxEZmhEuT8ORx/2OT/bzSeXFmGy6ps26qnqa4/Qe4E/A/hg9ek7hE3gSdCWJk6JMuvnplmmrbDp1M4XTPNDLRtyG0yvTXa6jiDvdMu1BLttfW14gOboOO3JpeZJjd9lyx3k5Zk1urzx3O6nqJkBVVUR+AUamGHfW8FUoXYBckEi6bWQiSWajRkkXqdr3MhF0E3x3JO+1hTy7GrYVSjpbG/gBWAms60hbH5jsUXZyinKTgdEiEgKw9mNd+Z51RWRDEdnH1ddAYGHbw88cPoGnQKYPaK4f/nRIPBMiz4V+O9stGdyemDFXntN8MB1dfqrraQudQd6+9O2N9v6OAFQ1BjwIHAtgrTGwKfC4iKwnIm+LiL3u9iRgoogUWrbjh1lpqOqXwFfARKvsROA7VZ3aVl2gEjhLRPKtMRwElACvZHpP0oGvQmkDHaVKSbdsumVIs1/3Q9BVb/B0yMf20mxvW2sjea9NyLEr/VXAJBH5DMNvE1V1sWUJMgHIB6Kq+pqIbIBZkhHgIYcdN8ChwEMicjxGm3dwy3hT1/0aY6HyvojYwTX3VtWVubvEOPxohGkgG5LLdYTDjmoz13W9kO7D6fTIjLnOs223u5B3e+t2ZFu5RDbRCPuJ6K5pln3Oj0aYAF8CTwOZSuGZ1slkApQOKJusbmfD9sgUx7kXcmkh0hPJe22DH8wqe/gEnia6C4lnU9ZGd5zwcD+4qcS3XBJ3JuUyLZvLup3RXnfA2nhNnYFeReDZkHB762dDth0lYfck/Xe2dTqCvNsDn5jaRo5jofQq9CoCzwU6msQzLd8eCTsZueQ6Nnq2cKpTctlPR1oX5aJeZ7fZHbC2XldHo0sI3LKrPA+oBTYAVqjqn6xgNNdilk0cBfxdVT/zqL8pcBrwE8aF9XxVTesl3l4pPFtkQ+J0Qp1U7XQVnP2ncvTxKp9J2x1RviPRncaSS/g68OzRVRL4RcD7qvoBgIhsbKUfDJSp6sXWyvN2KMiWVbpERIDHgN0t86BbgKOwIpClg65QpWRbL9s6Nrqj3tsLXg9wLom7M8rnqm5vhH+/skNXPd+HYTydzhaRa4DFVvq+wKcAlt1kI0ZCd2IMUKSqdp2PrXoZob0/mM78rE7HiaGtut3xAcl2bJnWybaPbOGrTjJDrsLJ9kZ0FYGPwoQJuB14H3jGSq/CqE9s1NA6klg6ZQAQkRPtqGU5GHMrtJdUO6ueu35XPgzt6b+ziLi7kXdvQI5c6XsdOk2FIiInAQdi9N41wOdW1kfATpaL61Kgj6NaGa0D0aRTBgBVfQB4wOq/lcdSLvTh7WmjPaoY2tGvux0vtFfFlEt0Jgl3R/Je218Kvg48e3SaBK6q96vqXqr6e+BtjCoETJSuHyw998vAdgCWDrwQ+M46H2WV/xFoEJFB1vkOVr2skIsfTlc99B0pSaf7SduRn7ntkdSz7S9b+OTdPvgqlOzQVZOYFwB/tqxJ1iO+5NEzwGYiciUwAjhSVaOWdcp7IrKrqs4VkcOBv4rIPCAPeKQ9g+kOkjg5qN+eNroTuoJIuyN59xb4Enj26FWxUNoit1yQX3vbyCUB9yQy76pJ5a6s25VtdxSyiYVSIqLrp1l2ih8LJQG+I48DXS2J2/XJwTicbeWqvVwiV+TUleTrk3fu0NuuN1fwCdyF7kDidhvkoB13e050Jql3N9M6n7y7D3xX+uzhE3gHIdeWIh1BtsmIojtZoHREH11d30cifB149vAJ3AO5VmPkUh3SGVJzd32YugPxdvS96a73vqPRW6+7vfAJvBOQKxK32yKH7XV3dAddeS7b6Mr2uzN687W3Bz6Bp0BHTCh2hE57bSTzXNuTd4c2urL97gxfhZI9fAJPA7mWoDtyybKeTObdwYOzI9vpqvZ7Avx7kB18Ak8TPUUN0tPIvLtZp3RkW13Rfk+Ab4WSPXwCzwC5lp47Qhp3t+9EdyD0nkSIPrl2Hvx7nR18As8QHUHi5LjNtvryQletRN8d++ussfukZeDrwLOHT+BZoCP12D1pncrugO6qN+8u/fQU+PcjO/gEniU6inC7msh7Crq77rw79NNT4Evg2cMn8Haio/TYPpF7oycTd2f31ZPgT2JmB5/Ac4COnIzsaVYlHYG1Je6IT97e8CXw7OETeI7QGRJzbyPztYW4u6K/ngb//mQHn8BzjI42DXT2Y2NtIfO1VQftk1Nq+BJ49vAJvAPQ2frr7mjvnS7WdhWGT0zpwb9P2cEn8A5EZ0njXv060Z0IvTeRqE9K6cO/V9nBJ/AORnewJulMB562+uts+OTd/eG70mePLiFwEdkCuBiYAmwD3KSqn1p5FwBlQD/gDVV9waP+KOByYA4wCjhPVWs7ZfBZojsQuRfWVqLpyutaW+9pR8HXgWePLlnUWEReBR5Q1WdF5EDgVFX9tYhsA1ypqvuISBD4HthSVVe76r8GXKGq/xORM4AqVb28jT7bXNS4M9GdxrI2wSfurkU2ixpn8mzG/EWNE9BVPLIEGGAdDwCmWsf7AZ8CqGoEQ+C7OCuKSD6wGzDZSvoY2LeDx5tz+A977uHf056LWJqbj0R0lQ78MuApERkPbAecaaVXYUjbRo2V5kQl0KDxTwevMgCIyInAidZpUwy+zcHY00ElsLytQjn6QabVV47QmX11dn9+X+3HulnWex0zznTQmb+/bo9OI3AROQk4EKgFxgCnqeqnIrIR8JaIDAKWAn0c1cqsNCeWA0UiIhaJe5UBQFUfAB6w+p/SWZ9efl89rz+/r9z0lU09Vd0r12PpLeg0FYqq3q+qe6nq74HhwCIraxFQYB2/jJHIbVXJesAH1vlwEclT1WbgXWArq84OVj0fPnz46FXoKhXKicB1IvI1sD5wjCVNfyYi74rItRgrlPNUtdqq8zRwIfARcDJwhYjsAYwAzu3sC/Dhw4ePrkaXELiqPgs8myTvpiTp2zuO5wLHZtjtAxmWbw/8vnpef35fPasvH3SRGaEPHz58+Gg/fHNkHz58+Oih8Anchw8fPnooenwsFMv1fhTGvHAccBxQBFwP/GilXaqqSzzqHg5sBkSBH1T1/jT6WxeYCDRgnIyuwpgxtunan06YAFf5IuBzq+z5IlII3AwssK7relWd5VFvd+B31rhUVf/cRj9jgb8AXwDDgBWqerWIVNBB9zHdsaZ7zSnavw2ox5ivbgKcraqLrbwOC9sgIn8CzlHVSkdam/cp3XtulQ0B51nXtgHm//YnEQkA1wJrrDH/XVU/86i/KXAa8BPGl+J8y4HOq69eF/6iR0BVe+wGDAJWAgHr/Hngj8Ak4GArbX/gnx51hwFfEp8HmAyMa6O/PIzJot3fYIwn6WvA1lbaGcA1HnW3AV6xjoPAbKBvG/3dAjwC3GydXwxcaB1vBHzoUacY85AUWOf/AX7VRj9bAb91nE8Htuio+5jJWNO55jb6+Ivj+CLgrkz+H+n8bz3q7Gr975Znep/SueeOspcDOzvON7b2hwL3WscVwCwgz1VXMI5tgxy/teNS9PUqcKB1fCDwZkffR39re+vpKpR6IIx5+wOUAt9hXOs/tdKSudrviYmrYM/ifgrs3UZ/W2F++GeIyCWYB6ya9Fz72wwT4ISIHGG19ZMjueW6VPUbYBMRKXNV3Q6Yp6pNbYynBao6WVWfdyQFgDo67j5mMtZ0rjkpVPUyx2kAI61CB4VtEJGBwCHAXa6sdO9TOvfcxmHAaBE5W0SuARa721DVlUAjRkJ3YgxQpNbXSBp99frwF90RPVqFoqo11ufb0yKyCPgFI9FVYT4fwbja9xORoCZ+HjrL2OU8XfIdGIkhnYmqulpEHgP6k55rfzphAgAQkfWB9VT1UhHZOI0x17Tzupx9Hwi8rqozRKSj7mMmddO55jYhIuXAHsBBjnZzFrbB6sNWXZwP9HVlp3uf0rnnNkZh1E63W6qoZzDSf3vuazJ0SvgLH5mhRxO4pcO7ANhcVSMicgtwBXGX/GqMdL7K4wFYCqzjOC/DkH8q1AAzNB4d8SNgJ9Jz7U8nTICNA4FGEbkY2BEIicjZabaRST8JEJHdMJLS2a62qsntfcxkrBlfjzNsg6r+XkT6AvcAx1oSabrtphW2wdFfEHPtJ2H0wUXW//A/pH+fUt5zV0iKGswcCVi/RRHJS/Pa2izTFeEvfGSIrtbhtGfDfIa95zg/H7iTJHpEzCf0COs4Gx14gj4RuA64lCT6PSAEDLWOtyWuK8zH6ArL07jGq0hDBw6MtvYZ68Ad9/J6jIpoCEbK6pD76OjTc6zWfS5r65rT7KMSeMzxfziorf8HJtSD/T/OSneLNbHuOE96nzBzKYXWcSY68CeAva3jdYFZ1nFSHTgwytpnqgNf5qhbiXnBSEffR39LvfVoRx5L2rgTo+OrBjbESI9NwA3APGAscLGqLrEk9n+q6kZW/cOBLTFWAbM0PSuUA4H/w/ygR2B+jAMxkv+PVtq5qlorIvsDp6jqPlbdCzCSWT/gVW3bCuUgjJVACCNBPoexyFiEkeauVdVZIjIAQw5jVbVRRH4N/N4aY7O2bYWyBfA+xsIAoMTq7wU66D46+m41VhG5EVipqtdbljitrjmD9r/ASMa25L1GVfe38jz/HyLyCeal8ZFlPdHqf9tGn+tgwj2cgnnJ36aqdcnuk4g8DTynqk9aViit7nmSfoYCfwZ+wMQNultNjPyA1W+9NeYHVfUzK/1HYFdVnWv9H8+w+qogtRXKgcDBgB3+4r9qPKo77D76aBs9msB9+PDhozejp1uh+PDhw0evhU/gPnz48NFD4RO4Dx8+fPRQ+ATuw4cPHz0UPoH78OHDRw+FT+A+fPjw0UPhE7gPHz589FD4BO6jV8Fy/uo27fjw0R74BN6LICJbi8h7IvKJiFwlIneIyL+smNvJ6vxRRFal0fbZrvMp7SW5dMcrImUiUmt5lCZrq1hEbsd4HOYCw0TkRhHp0fGEfPRs+ATei6Cq/wPeAz5R1atU9SyMm/6BKeo8DqxOlu/A2a7zrVQ1muVQ7b7THe8fMbHgT0zR3J0Y9+9lANaLYYJ13F9Evs1wbPOAzzCxanz46BL40oOP/liR4UTkasxvIoqJGXKjs6CIlAJPAx9ggic9oapvicjBQLmIXAXMwMTguNMKcXqtVfYIIIJZoOKvqvpcW/21NV4HxgHnANNFpNVKL9a49wVOcCSvgwnyBLAx8E2K61sXE1TrG8yKOtdY8VheBe4RkStUNZbG2H34yCl8Au+d2FpELscsKPCAqr4tInsC26rqHtAiob6hql866sUwgZnesoIuvQ68parPiMiNqnqVXVBEzsUQ9gnAVxiyDAAfWOSdTn9Jx+voZxtMhMKlIvIaZpGDB1z1x2KCY6lVZySwwEG6G2OCNHlen9VvGLgbGIoJnoaqNlhqooGYYFs+fHQqfALvnfifql4jIh8AN4jIIxgSK7biVwP8THwFFhsC7Coi2wHNHvmtoGbhi1cxIU4LMKFdSbO/pOO1yRiz+k2diGyCiUJ5Eq0JvADzMrGxCYawbWyBkbyTXd+DGAn8Q2AmcK6jbjNmDVYfPjodvg68F0NV38eE4T0QIyUvVdXrVfV64GEMWTlxPDBEVa8BbnXlRcVgE4+u7saExd1UVadZaen0l2q8WAs1VKvq5ZaO/DiMKsc9mfkziZOXmwKFVhvjgN9i1CPJrm8bzGLK22CWFjvSqiuYZfwWphq3Dx8dBV8C70UQkS2BnTEr/PxaVd/ETMI9Yu2/EJHrMEtt9QMuFpE/An1F5GSMSuH3InITJr52XxE5SFX/g1ns+Warn3cwy8+djIln/Z2IrAHetceiqm9YViYJ/aU7XhHpg1lnMuwoP8Ea100icqqqzrD6WiQiC0RkqKouwEjgjSLyFUYSnw4cBfzL6/owkvmtIvIjRiq/1+pyE8yCIo1Z/Dt8+Gg3/HjgPnoFrC+D0zEvlRmYZfjWpK6Vsr1CjGrlclWdm5NB+vCRIXwC99FrICKDMV+db6vq+Ha2NRSzUO/KNgv78NFB8Anchw8fPnoo/ElMHz58+Oih8Anchw8fPnoofAL34cOHjx4Kn8B9+PDho4fCJ3AfPnz46KHwCdyHDx8+eih8Avfhw4ePHor/B8pHmAcUldhdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot polarization ticks\n",
    "im.display(plotp=True,pcut=0.1); #pcut controls the location where ticks appear\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8f727a06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAEZCAYAAADv+lzvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAAsTAAALEwEAmpwYAABgjElEQVR4nO2dd7gkVZn/P2+Hm+ZOzjDknDMCkkUEFUUMu4IoiqvrGnddw66KYXVlVdxdw4qYCAr7E4yLKEmCItkBZoABZmAIw+S5E27qVO/vj6rqe7q6qqu6b/e9fe893+epp6tOnXPq1OmuT7/vSSWqipWVlZUVpMa7AFZWVlbtIgtEKysrK08WiFZWVlaeLBCtrKysPFkgWllZWXmyQLSysrLyZIFoZTUFJSKLROSHIvJgwvh/IyKrROT1gfAzROR/ROQLIvL51pR27JQZ7wJYWVmNi04EfgMcHhdRRPYANgAvBsJ7gMuBg1Q1JyK/EJFXqertLSjvmMhaiFZWU1CqegOwwwwTkYNE5GoR+YSI/EhE9vTiPqeqd4RkczzwvKrmvON7gNe1tOAtlgWilZWVrx8Cl6vq14FrgMti4i+gEqrbvbAJK+syW1lZ+ToUOFNETga6gf6Y+BuA6cbxDC9swsoC0crKytejwC9V9TER6QTeFBP/XmA3Een03OZXAv/T6kK2UhaIVlZTUCJyCnAhsFhEPovrHl8MfFxEVgKLgeu9uAJ8BtgN+BsRKajqzao6KCIfAL4lIhuBxyZyhwqA2NVurKysrFzZThUrKysrT1PBZVaAdrKE26ks7STXM5u61x8HNXTDIpL0B3yzqp4VkUcW+CBwrleODuCzYS63iJyKO95xnRH8HW/oUFM1FYBoZWXVZCX581DVeTVO7wx8FDhcVbeJyKuB34jIfqq6JiT+pap6ZUOFrUOT3mVW1bayyNqpLO2m8a6b8b7+RJKIxG4x2gFcoqrbAFT1VmAYOKHFRa8payGOoewDFy9VnYqu64RTKhVvSzmOE3lOVTfjDv4Gyj3ZHcDGiCTnisi7cJl1M67FWKyjyIlkgThGsjBMLr+uxgOMFsjxSmgBAswTkYeM4ytU9YqIuKcAzwN3h5zbhjvm8TKgC7gRmA18PHGhE2rSD7vRNrjBNijChNV4wWmKQLGhm0ylUprNZmPj5fP5h1X16NhCiHQBfwQ+rKoPJ4h/NnAD0Nvs59taiC2UBeHoNV7WorUUa6tZdeO5ylcA/5kEhp5eAHqA+TR5quCk71QZL403DP3OpGZv430/Y31Nq3ClUqnYLaEuAx5Q1etFpFNEdg1GEJGPeFakr4VAHtg8+juplAVikzXWD+5Yg2u8QTke17OqVJIe5iQWpIh8CtdLvVJEeoG9gPeIyDwR+ZOIzPGiHgm8zUuTBj4MXKuqpWbfmwViEzVWD087WGxBjTUg2+3+p5pGC0QR2Re4FBduO7ztce90N7A/rlsM8APgb0XkDuA+3J7ojzb/rmynSrOuMaHzHwu1sj1urNr6JmGbYkM3lE6nddq0abHxduzYkahTpZ1kO1VGqVa7p5NF5r00GyzjOUxnqmqy1rUF4ijUCmBNJghGqVVwbDUYLXhd1TEOccLJArEBWRA2T62Aox0y03pN1vq1QKxTzQbXROoxbfVD0EwLrJXWnAVusql7E1EWiHWomfCaiG2PUfm2c5tgq+A11aE4We/dAjGhmgWZiWxh1lOGZlp5o82vVdbiVIWibUOc4mo3GLYDBOMULONoH6BmwKcVAJuqULQu8xRVM+DTLnmMp5ph7TXD0rNQbI4m6/1aINZQO4BsooMwTKOF42jBOBUB1mxN1vqzQIzQeMNwIrQ1jnc74WjA2GwoTiXI2jbEKabxhOF4g7gZ1xlrl7hRGFkoNi7bhjhFNF4u7mRyrUfbodIIGEcD00bS1cpvKkBxst6jBaKhiQbDdoJgLTXqFjcKxsn6sLaTJmsdWyB6migwbIfB4c0YE1hPPvWCsREoNhOkkx3KImJd5sms8YDhWIJwrDpoWu2u1hPfWpet1WStpykPxHaH4VjAtllqdMZKvVZjPeCqF3LNguJkh+tkvbcpDcTJBMN2bU9sBHZJ4tZrLVooNleT9b6mLBDHGobtBMJmw7NelzdJmnrA2M5QnIyybYiTTJMBhu1kPTbiKieFYxIwtjMUJytYJ+M9wRQEYrvCsJnx2sF9rmcsYjOg12yLstH4U0WTtU6mDBDHAxJjDbnxGMLTTHc5DmrNshbHGnKTDart7jKLyAHAe4EDcd/g9wLwC1X9TVzaSQ/E8Rq31ywYjhcsk6qRITijBV8Sa7GZUJxsQGuG2rU+RORtwFuBm4HbgAIwBzhVRF6nqu+rlX7SA7FZakcYNgOmcWoUBvVYhFFxaoFxtNBMGqeRuK1I325qx3sRkRSAqr415PTPReRQETlIVR8POQ9YICZSu8Gw0XONxKsnbhLwRcWLg18j4GtH93myqB3rTFUdEXkh7JyIvFdVfxiXR/s2BLSJWuFqxgEt6nyj58zzcfFGo+A1mnkfcWlqlalZGusOuXaV34YYt42TPhYMEJHFwD8nSWyB2ESNti2vkXNJ4FIPKBvdkuTfyLmw8Kjr1CpDLU0WUI2l/DURa23jpNki8iGjnPOBW4FvJElsgVhDzXYrWwHDqPAkkGy29RSXd7MA2GwrMsn5euO1Kn27qF2BqKqvAfpF5HUiMgsXhj9I4i6DBWKk2gGGjcCjURC1wkqsdb6ee2iWtThZYNQOaleXWURer6pXAkcBdwD/T1X/W0TOTpLedqqMUq2EYdLwZrmWSeMl7TwJxvXPx4X54Y2GNRIed66ReJNVzbAARSQLfBA4FxCgA/isqt4eEf8C4J8ABe4C/lnDf6T/ISLvBNLAbsARIvJzYH/g93HlskAMUTPdp2bBsNlhcWWrpSTp/AfGjNtoWBCg7QDF0WgyALUJ5d8Z+ChwuKpuE5FXA78Rkf1UdU3gWgcDlwEHA1uA24F/AL4bku/TwO+8/d8a4YlMVgvEgJrdrlbvuUYhN15wDCoMcH54o9ALO45L02pNBqiNRk249x3AJaq6DUBVbxWRYeAE4PpA3IuBm1R1k3ftHwOfJByIn1TVZ0LK++ckhbJAbFCjaaQfDQzrPU4ap1Z4UgWhFxVez/FowBkVVis87pyVq9G2EarqZuAa/1jcCu8ANoZEP4ZKa+8J4CAR6VbVISOPFDAr4nqrRGQuMFtVV0aVywLRUDNd5XrTNhuGjYCzVnirFAW/KBAmsR6bAcUkGk36iQzdOtoQ54nIQ8bxFap6RUTcU4DngbtDzi0EthnHW3HbHecBL/qB6g7MPk1EXoML0DVAEXfq3vHAG4B31iqwBaKnsXCVmwXDpOeaYV3Wo6gOknrVbEhGhdXSRAbWWChh3WxS1aMT5NUF/Dtwkao6EdHCfpxVhVDVr3kdMN8FDsK1Ol8EfgW8V1XztcpigVinRuMqJ43fCPBGsx9VjlrhQdVq40sClzB3OwjGOEg2AsXxtBInspo1rMZzla8A/lNVH46ItoFKV3gWLiDD3GtU9WfAzxopjwUizXOV6203jAurF4xJw5Lux5XVVBgUzDB/Pwxy5mdU3mFp/HLVA8V61GrYTWSYNrHclwEPqOr1ItIJLFTV4HzkB4H9jOMDgcfN9sNmacoDcSzazJoJw6RQjPr099PpNKpKqVSqme+mmYDC3G2166mvN0OHQu9gqRwWBUT/s14gBAEalk8z2xitwtWMcYhePp/CZdCVItIL7Aq8TUS+g+vivlFVtwA/BG71OkX6gIuAy0ddgBBNeSAmVaPW4WhhOFoAmp/m/j2z17FDcpy5aZeqc+Y17shsJI3wRmdeZNlFhGsz61hMJ+eWppXDzPMAmUyGX8zYwr7ay8Hbs+VzJiTrAWUtQJrXHQ0Uk5RlKgJ1tC6ziOwLXOodftg49UXcRV33B3qALaq6XET+GbgFcHA7Xv5nVAWI0JQGYrOsw9Hm02wAJgkbIM8QBUqlUk1r0kHJkqqIF1QqlaKEkkLKFidQBTrHcXhKtzOTNAeUUuXwYDwTdHFudJyFaNZfM93pZqldylGvRltmVX2akE4RQ/MD8etuFxSRc4AZwCPAc6o6GJdmSgMxqRoFXhLrcLQwrLVf69gxABaVh4h4oINisRhZD+l0GgdFvHi+wmAHbmt4qVSqAGDYZkIx6QMYZTWa9xSl0ViJU03tXh8i8jVcqOZxxy1+FXdmTE1NWSAmhVwrXeVGYRgHwbhwVaWEQ9qw/MLSiIgLThUKhULN+1AgpZTjhcEwlUq5JoEmA2LUFlSYtTjWru5Ugmaz2hBbrK2q+kkR+ZSqLhV3bGKspiQQm+Uqj+aaSUDaCASTbo4oaU2VLb+wrQxEXCD613WcyqFijuN48SCfz4e6wD4QwbUQi8Vi1bmw/WbDMc51tlZiMrXzS6Y8zfM+/QdqflREU1MSiEnVTOswST5JYRgW5jhONfSMsOB5QciQKoMumA7cH72DIuqCzryuKd+19uNFWYepVIol0sNMJ0OhMFwBP/PTcZyqcBOocaoXiuOtditPEk2A8j4tIk8Ajoi8Bfh+kkRTDojNsg6b7SoHw6Jc5UiLzwNe0s9Tnp+HqtKf7w8Fpg+lE4bm0lkShoe3h5YVXPf33P4FdBeUXC4HVFqHPtQKhQJnPJUB2U7OcYf++KA0oWlC0YdFcD+JteiXIw6KSa3EOE1EsDWqCXCf1wJ/xu2tXqaqTyVJNKWAWA8MWwXOMBjWOh8FxDDrz/w0NzPMbzMcHh6uCA9LD6B/HWBAHdBweAMMDw8jj/WzHUiFuMlBKEpKSKfSZehFbeb5IBjrsRwbgWLUdzkBQNBy1WOtj6OeBs5Q1RvqSTRlgNgswNXKqx4Xu5aVGOUa1wNBE4A+BP1jx3FwZu3GsJOme9sqcrlcVT7+fnrBUTy5uZdjpz/Kpk2bgOo2xFQqxcyZM0mn0+zYsSPU5a2wBNNp0qlUhYXob2Fh5ubD0IQihC9pH6bx6HBJookG2wlQ1l+o6nL/QESOVNW/xiWa9EBsBIStdqvD4oTB0FeYOxsGwlKpRDabJZVKMTg4SLFYrICg/1ksFrlxdQ+PbhA+u1ue4eHhUJiqKqvXZfnFo8pBxw4wNDQU2Ya4MXM8f1g2jQ8d9xCbNm0MdYNnzJjB87ofi3v6KW1fXYaf+enD1ASjP6vG3HwLxQRiLavFB6G530jb4kQDV6s0AeqgU0QuBZ7E7Vg5B/cF9jU16YHYbDXTOkxiGUZZhSYEzf3rBubzy3UO1+02TD6fL1uIJhxLpRKFEnSmYWBggOHh4Yp8TCDukDQzu4WtW7cyMDBQcR8mjDaQZSinbN3ax9DQUKjVVyqV+M97urjgyDQHSr4Khv5WnLaIX66ZzTv22Eomt821VAOWow9JH2omIKPgWC8UG4HfVAHmBLjHo4BfA7t7x3OSJLJADKjV1mFUm2LQKoRqGNZyiX3obS1CwXHb9fL5PMVisXzOjJtPCx1p6O/vL8dT1Qr32nEcdmTTzO6Bvr6+cocJVM9C6dMsC2YJfX19lEqlsoWYyWTKAMvn82TTbvnypXzZIvQ/S6USmUyGggPXPulwysIOdtFiBfRMazHf1cv5zwlf2TXD4Znhchy//TEMjEHYhUHRqrYmSBviB1T1Pv9ARPZJksgCsQ7VC7la6eI6T+JgaMLNBN5ACXozwtDQEMPDwxSLxaqtVCqR74WOtLBj+w4KhUKVW+23E+7oTjGrR9ny8paK8pttdV1dXWzaIcyboezYsaPCivPhmE6nyeVyZNNCsSTla5qWoQ+7WewgJTNZM5RlUbpQDvc//fw6tMSyfmFNPs0hUqp4SINDeIKqBcE4K9HCs/0tRBOGnk4Dql4tEJQF4ig1Wlc5DIxh7YVRICwUCuV2wWKxyIBWA7FQKJTP+yAqdAudGdixY0dFfuZ1enp66BtU9phTYnBwZBposAd5xowZrF+n7L0wR25zrgw3Py8fZL6FWHSEYrFY7hwJbrnBfpZMF9YMZjisu1gBQ8dxyGQyqCpdTpGeVJbNpZEpiME2xjhrMcyFHi+N9/XrUbuXU0T6cFfGEdxB2Vtx112sKQtEQ7XgltSVjrMig5/+fhCCUW2FpkVobj70BtIuEPv7+ikUChQKhbJL7MdzHIf8LOjt0AogmtdRVbq7u9k6qPQsyLNjYKD8EPhg8dv1Oju72LAderND9OfzZDKZ0PbIQqFARwYKHhDT6XS5DkyLuL+/nz1mpXixH4rZYgXk/DSpVIpcLscps3uYmabm4hO1rMUkLnS9oJpIYGtEE8Rlfp+qXg/lFbkvSJLIArGFSmI9Rm1RbrLZJugDzt9KpRLTOoVpKRgaGirDsFAokMvlyhaliJAvKRlRBgYGKiBrliGdTrN1ELpkmA1DQxVWoe8Op9Np6JhNsQTdsp1txmo3vkTclW4KhQIdnoXoW22mtemHFQoFztolz3BRkeKI9efD0N8fGhri0p48kk9RyqRD/2yAChg24kJHfXeTGXpxavd792Ho7Q+LyB5J0lkgemqkMyVJW6EZXqsDpVbbYS2r0Aee/3lB8S84jsP6wUHy+Ty5XI5isUg+ny/HyWaz5EuQEafKQvQ7VgAk3UF/DjKO2xPtg8ts83Mch2zXTPZaJHTQV160wZxp4t9LoVDgwkM20yU5+p6rhL+fxn/QFm1+0LUCjeE4fj2ZYHRd8gxQGR6sZ9O6rWVJ+ulGayVOdrV7XYjIHYzMY/aXAIuVBWICNQOWYeFJrMNgp4kPQxOE5rZ169YyCP2wXC5XEbejo4ML5q9AC8Pc3N9f5dr68Ck5yt+fVGLW4HpyObdt0O85Ngdnr1h6GyfPf5T1q3JV9232RhcKBbY/czv9qRSZTPhPL6o+fOsuvC6LQLoqL7MMPgj9PIIdMGZ+tvc5XhOgTu5jZFXtHequvB0rC8QG1ah16O8ndZeDHShmu6AJvODmD7vxwehDsVAo0NfXx41X/XdVB03wHv581+2k/3wnkkqRMYbIgLegg9eDPDAw4PYgZ7Nks9UrYZvjB4OzV6JmlYTVSb3teGHHJhjLq+8ErMgoCAbDa8GyUZBOBABPkDbEX6vq8wAicqiInKuqP45LZIFIczpTkqQJCw8Doe+2mh0dwTZDH4bDw8MVEFy86CByuSVs234nzz77dJV1aPY4B11l0wJLpVKI366XSoEBjUi3MwSAYQOq46bnmZBUrXRhg1ZiLQXBYh6b9xBl1VorMVoToD5eA9wPoKqPiftq0lhZIDagpJBs1DoMgrBQKFCig1KqE9XNVZagufXv2IPrr1vEuW/tZHBwsBzPh6EPxIvO/jH5wS5+fs876e/vr7CSfDB1d3dzyjFvxyHPPQ//KhSEIsL06dPZc+c3ke0a4NkX7qhoY8xkMhX7XV1dZcj4YDSvWcty9K9dDxSDZfU/zfZF857qtRKnqtq1DkTkXbgvodpNRE71g4HhJOmnPBAbsQBr5ZHUOqzlJocNs3lq6yK+/n+dXP6OAnlv1ogJRH/c4dBgN73TYfOWtQwODpZdZ3/ojQ/EDatmsmNThoGBgfKaiMHxhTNnzmT4sbcxf/cc8KsqcPjxuru7+esdB7L/YcOk03dXTcfzZ6vMnj2bu+86k4MPydHbe2tFj7UJw66uLja8vDdz5g2g8kJNtzSJ1eiQ5tltPew6q0hPthRpKQbHJfrXSDqIe6qozV3mXwN3Au9jZNxhCVibJHHb3lU7KMrFrSdtlHVoHkeNOzTd5PXbM8ydDrnBrWUYDg0NlbfBwUF3G+hizhxYu3Yt/f39DAwMlM8NDAyUwekUU2Q63Kl7ZvuiD0zHcRgaGiLTAU5xpN3QhKY5lS6TBcephKBvHaZSKbLZLOl0ms2bYWhQQi1IP253dze/+fl81r60GJFqNxpcK3HatOnc9j9HMbBm78g/FMdxGNYePnxDN89s6anqQIrbrMJl/haitvGQqm5T1edV9TPe5/Oq+hJwWJL0U95CHK2SWodBCMZ1pAShuG5bhl3mUe5F9tsOh4eHy8DL5XL093cydy489ezLZevRBF15kYdiikyHMjg4WP4B+0CCkYHP2U4oFaoXTqgCYgackpDurJyKZ1qIrjUGCFXtiFVjEhVEzEUbqt001RLzlxRJd+Qr6tEvq78/vTtHOjWNDf0pnLkj4x/NuOb3F+aqJ+1MqZVuMqnd70tElgAfwX2VgACHAEfHpZvSQGx2Z0owba2e6Frth36Psg8wgFxBWDLXKa9O47vJPgyHhobIZDJs2QKLF+fLlp9p/fkwTKfTFPPQ0eNUvOzJBVumvF8sFsl2glOqXIy1ctku13rMZKFUSlVZe+aniHDAASXmz69c6SbYbpjL5Xj3B1bT3TNMyal+8HxIDgwMcMjrl3rlHel0CY5pHB4a4Jq3O/R2lHCc6jf7mWMmTbfZl21XrNYEuN9LgV8AZwE3AG9LkmhKA7FeNaPHOc46NDdzQPWZO98JwIYN4R0qQ0NDLFiwgDWrYfc9XECaA7PNcYwdHR3M2X8FknZBaYIpWN6Xer9EKqWkNlR2epgWXrFY5OAT/0QqU2T9xuoeY39/cHCQ/Q+4jZRxveAwHHBfQNXVsx5HqZjel7Sug5ZfsVikMz2IU0ojpCtmxZgdLKYlP54PfLsDts3bEH09oqq/EpF9VfVWETkuSaIpC8RWdabU03YY1ZkSnJGSzxe8tr+hMgjNtkP/eN269Rx+9M/ZPtBXnrpnTtnzr1coFPjy995Q/mH7VmKw00REeOSx++ns7KCzszPUqkun3Qdj9UsP0tnZQUdHRxUIo+AXPA8Y+1J2k3021NN+a1qL5r5pPQbj+FvYH4NfJ3GdLFNFE+DejxKR3YD5IvIO4DTg3+ISTVkgjoWSWIdRbYjVYw+rB2EHe5j7tvTx3HM/Lrcvmq5ycPC12fDttwMGlUqlQKiCWnDzwZi0kT343hWzHOZx8KELtvVF1bkJuDBLMDiPOuz7mAAP/LhqAtTPfwG9wPeArwPfSZLIAjFESXqXk3Sm+OdqdapEzUoJwjCXyzM8HD72sAzF3HDFQG1/xWy/7TDMAowaZO2fz6RGepijBlyHQS8MZrVmpox0npigglQquoOjVidWsK7NttKKzps6XOdaoJxq7YsTwGX+JPAVVX0cOC9pIgvEJiqqM8U8H4RhLeswOFXPBeMIFE04+vulUqkMQ7Pt0ByrZ8IpOB2vlnVXy+oz806qJBANC6v1hxWMH3SDHUdJpapd6rA/rrBmhGDekxF2cYr67ttMjhovlRKRtKpWL8UUUNtjvhVqRvthPdcwH7RgWC3r0IdiLpeno6OT+bMPRSQVOUsll8tx4Rs/x9+94iFec/p5VUuDmduXP3ITX3jX/SxevLiq7OaA6Te+6uO89uQPVllw5n42m2WnhQeycP5eieojTGHQTaWiXfSgxRnmase12Yb9OQWh2Mi9TAXFNY8kBaaIHCMiK0XkohpxThWRFSJyp7G9JSbre0Rkf+P4U0nKMyWBWEv1/tiTuMvmfnArlDoYymdjLcTZvQfy2++dyN47v7HCQgyubJMd3psVf8mwrX992WUOg8Dw2t15+bH55WE9QSABZDIZtj1+MrnVJ9e05rq7u1l626vZvPq0mnXV3d3N0qWns2NHeIefn2c2m+Wl9buTKyysKFMtKNaa7mdCD8Lnjod9V2Hfr4Whq2YAUUTeBPwjsC3BJS9V1VONLe59y/8G3CQiz4rIc8DHE1xj6gGxkR90kocizqrwzwUf0AefmcVH/mch6WxPpIWYz+fpW7MXnd2wqf8vFSD0931AMrCI+bvCk08tr8rLbEcs5TrongHbt2+vKKMJRrezxR0gXcu1VVV36Cvx7uyqlSm2bMnUrKvOzk5+9otZPPvC3JoP3PTpM/jeFw5nw+r9az6U6gjP3LI7+c0LIju0ov6wgt/bWKndwVtrcQ5/S6AHVfV8YEcLivhFVd3T2/YA3p8k0ZQD4lgozF0OPlz+QzicV044WCjk+qvWPfStw66uTp5a2sshxzg89fSyKhj67YWZTIaBjdNYsLvD888/X9F+aF4zk8mQH8zQNc2peJOeLxOKjlNexKbqh24+syKAcTyaB1pEKDlgPlNhsCuVlIOPKjB99nBNC6Wzs4O5uw+TnbO+7il7cW5zkriTTUmswyQWorpT6pLqXBG5Q0T+JCKfFZGa/R+q+g0RmSciB4tIbwKLErCdKi1TrXZDE4on7r+G4/YuMjxciFzqa+7sRaTSyuK9XmLFAwMV6yCaq9jsueeerH9O2OsEd0WcqDLMmDGDoe3QPTNfDg9zSRFBHSBVu7NEVVm0BGbNVbY5leHm9YvFIqefvo2584Yi8xIRBgcH+exHV5FOO6gTPQwnlxvi1HNXeNca6Sn2r+2nyRdyzN5nLf7/fxIY+vt+nqPpRBht+nZUwvuZJyIPGcdXqGrsi55CtA24F7gM6AJuBGZTww0Wd+zhZ4BlwLUicpCqfiXuQlMKiHH/4En+4WtZBLWsiSgXzJ+REtax4g+ofmbVCli0lOXPlSrWNgyunD00OMSC4+9lna6qsAjNH6/jOG7b4p4/Z1M6H1peX8VCga2LL2GH0Qvty+/0AHcKXefsb7OlkK05s2R4eJjZcx7AXd06Op47lW4Ap1QNw8oyAAiOo5jPZxSAHEdJp0c6XJJuQShORsDVq4Qu8SZVjZ0/HCdVXQos9Q77ReQ/gBtE5J81+qE9TFUPEJFPqeqvReTgJNeaUkBsRPVCMu5cVPtV2IBs10p0LcXh4cHyS+WDUPR7kVc8tYLHlp1TzssEsDkbZPv27Xznmn8km82W5y4H5brLDitWPE5nZyedXZ3l8DA5zsgrQv144Z0cEMeSIHjDOnwq41cv3RX2nYi413cBWtsyDJOZ/1TXONfBC0AP7utFN0TE8Ttq/C8z0RzQKdOGOFbtO3GN8rVgGAbFfL5AsVgogzG4ck1waE1YXlHtY0FFuaZxiotf3SETZu2Fjzmsdd6/lbh2rGqr0swj3OJPWmdTUc1qQ6zjeh8R91WivhYCeWBzjWQLReRy4JUi8g3cWSuxmjJAHC/VeshUterl8OY0OxdyldALQtF8DYDZiWKCMDhLJdx6Ch8ULRFjASvaGqFq38wjeD4YJyo8WBZT5nJiUWAMuy/3/h2S/GFFfY/BMNux0lwgep0hfxKROV7QkXir1YhIGvgwcK3WHmj9MeCvwIvAU9hxiK5a9Q+ftP0wWIYwSJlWnQ/JYFtiqVQ5FKfi9QIGFIMwNK+Z5MEN/THLyLlg3GB42AMRBKTf9hj10MS5yCZA02l3HnUUjOMe0loucxjwrFw1Y9iNiBwlIncChwOfFpFfeqe6gf1x3WKAHwB/K+6rRe8DNgIfrZW3B8s/466efY+qJnrvhG1D9NRoW2G97Yf+pw9FcxhItcvs77vW4EnHvJkHH72ZbduerFjf0JyvfOrJZ/CGva5mRe5yvnvV54Dq+bgA//TebzJ921k8vOPz/Pm+W4HqV3OKCG887Z9xhmexbO0PGRgYKFt0wZkruy05GpECm/pWRkLJDRvJP3gtf9PAHGYzrtkWWsk3rZibPJrvJ/hd+fnZjpURNePeVfVh4NSQ8Bdx2wf943uA19aTt4h8CngnsArYW0SuUdWvxqWb9BbiaFSvZRD1QPmfYVs6laVYyNReMbvksMdu+7H69+dy/E7/XV6yK/gKUcdxOOWQd7Ps/2YwvWdOhbUTdJvn8ApW3DaLEoNV92ECLffyEbxw//7l47Atk8nwxF+O57lHT6pwpYMw7O3t5YGHT2fb9lfUjNfR0ckTq3dnsLCwBgwFkRS5fA/QVeXCh5U36juLc5fD4k9lJXGX2+DP4iBVPUhV36CqB+KumB0rC8QmK679KfjgrVy6B9d94xA6O2ZGrnyj6rBr5u/oWw8DM37N0NBQRceLCdB5xdOYtQhuvveHkQ+4iOBsXcLO+8Pyxx8rhwXbAbPZLANbOpi9UNm2bVto+6H/UqjBAejsql5RJgi+51YLfVszkQ+NiNDd3cl1N/by3JqZ5Thh7Yxd3b38y7/vzMrnF1Rdy4ybTU3j3iv3ZWjDoprf3VCxl7//r7k8u7431JUejRpJ387gbdJMlVZqReD4UQARObBWonEv9VRRGJhKJYfeuVs46bUFcvltFe2AxeKIG3zMIWdzz/UzOP7NO7jxth9V9R77aQ484EBeeHAGexw3wEMPP1Sxwo2p3XffnbVPZZmzxwB9fX0V50wLbN68eWxZK0ybM1QBtaCl1tXVxWA/dHaXaq6b6G54s1+i2/pKpRLTeyGXT1WVzUyHFjnhGIfeaaWq/MwH0iHPQWdsomv2yJTZsD+tjOR4/QlFZk0rVJ2zqtQEsBBfKSJXi8gXRORq4AgRuQT4Vq1Etg2xhYrqzDDDZi1cS/es1QwNFSvg5i9EUCoVcXLd7H2Ew9MD/1Ve4iu4YEOpVOI1J13Ipt9A73F/Cl3aytdhBx3Lhrth8fGrgeqOEP9z3rz5bFsG2d7+mlZfNptllz1LzF+coy83YtEFwVkoFHj92RuZNStXjhP2EOVyOT7zvpUIDo5TGc8sXy6X423nrEbEH99YOYDab1N0tMjMXSrBH/ZdpCXP647Z6g0uT/7qgqmoNgBenPLA7d7+c0b41lqJLBATqlYPbSMdMu6DS1XnSvW7VZR7HvkZqfR1vPjoi1XxzLUOf/qrb3L2SZv4+R131izTQ0vv5dRXfY17Vr5Qdc50hVetWsmep17Cyh2VrwWASpdp8+bNTFvyYzYNZsuvFQ261X6754zpT5KSNCLRliSAUypEutTmpz/sxoShee8+FGt/N2Fh4U0fwTJNxc6VNrEA4/QRr3MGABGZr6obReT2WoksEGmNSxTVI+040UM8RkBnDs1xKDkltvRtqQBgcHiNqrJ27Vp+cN2/x97PS2te4H//77/p6Bh5BwpUu0H5fJ7nnltFV2cXnV2dFatlx7UVhr1HRUTKrxqIW8vQtPLi4ni1W2UhRn0PQfl/TP5n3PcYpqkGxjZoI4xTn7jLi033js8B3qoxC0pYIEao0WE4teNXpg3rgQYMK7Ayjm/pRPWMBnuS/bzCfrymOwvhbULlxnEhEnpRm3neTOsvwpC0/SkImXAYRluI/r3GWYl+vZqfwXNTCXhxmgB1cSPwCOC3lcyJjjoiC8QxlMiINeIrDIxQ+RkFwjhFPcRRIAwe+1ALe39y8EXzQYswCMLKdsDk1mFUGc06Ne/Fr5daVmLwD2Kk/uPr0YLR1QSog5Wq+jH/QET2SZLIArFJihu75n9GWSEjVp56D2b1sA9zELepqPAwmbDp7e2lUChEwlDE7TAJAikKdlEgrGx/FNLpeOswzFU2yx/c9+vIBFe9atYzPtmhOfJdtrVuFpF34w7MBrgQ+Lu4RG1/VxNRYR0wUW5YRA6heey5xz6cecr5FeHBPLq7u7nu39fyL3///YqluIIP6HGveCXvO+4R/uEd3wAqF381QXbKcRdy4MzvcNhBJ5FOp0Mtw1QqxYIFp9Az7aTy+5vNtkYz3rRp3XR0dCRyl02FWbH1aAI8wBNK9TR7jJPeDbzB+3w3cESSRPZX0kIlsVIymQxZmVc+jvodLV68M4dP/2+GHvoIJ55wemR+773gU9z7w0XsPffVVe2HJlBOPvh9PHF3ilyhekhNRfz+g3j8/iwFZ1sFNINAXP7Erjzw0JKaFuO0adP42u8O4aGXDgx1s/1rp7PT+MOTu7G1ML/mw6WkWLN5JgWnO7aezfuvJTvkMJkmABA3qeqbVPXdqvpuEliHYIE47lrxp4P52Zf3J+WtTuQPITF/UPPmzefwGV/n0T92cMS5K7jvgbuB6ngdHR0cPu3j9MyEXz48sphwEHIdHR1kN7+SxXvDH++/KtLF7ezsZMPquexzqPLCi8+W4wStxO7ubjZsFBbMd8oQDrMORYRiEfLFqLZFF7RdnWl+dV+GZ9b3RLYvAqQy07nkR7N45uXZ5XoLU5rp/Pjz+7H5xV1qfhfFfCe//tECNq9NtFLUlNYEAOKjInKaiOwqIrvi9jLHygJxnLVo73Wcc/EgmY5CuQfWVE9PD4fN+iIP/q6L0y5+nu/f8H5KpcpVj3yoXPS3H2XZTb3sd/bL/O4Pv4q85umnnMUTd2fZ5eh1rH5+dag7mkql2G/fA3jmMWHhnlvJ5/NVgPO3mTNnsm49zJ6dr7AOg+AESKfB0cqB26bVKSI4pQK7zBf6BipX6g6WMSODnHEMZDO1X7ebyuY4/a3bmT5vU8146WyBV7+1jzmL+mvGm+qK+h0Et3HWPwOXAFd527uSJLKdKmMs3yXzH+re+evIzXiewcE8In5P9Mi/7PDwMOz8NGf83U5ce/tHy2AK0/H7voMVy+B3j3+66pzpkp6w/8UsfQCeH/5FlfVobotnv4InByA7cyWptdWA87fOzl7mzYW5cwcRrY7nP0CO4/DWVw4wb0ahArxm+VwrssgXz3sZtAhEd6oUi0XOP/2FinoNU6GQZ5cD1kQ2R4yA1qF3ZgF3yb3acae6JkA9/IuqXukfiMgZSRJZILZQYZZX8Idkng++mN3vlb1z6bcZGh5ieHi44lwQYO//0jG8/uw38cvf/L/Ia6ZSKX7wqw/x+je+n1v+fF1VWc1t6Ypfc/yFm3hhw7M1xx6uX/8Sp5z4GzKpLhxHqjpf/H3Hcdh77movL9NqTFfVFVqoqp84NWOAfVgW4/nwt2uPdTuWyZSqXikih+IuI/YUI9P4asoCsQUSqRwcXGnhhLuAwc8K4KXcBQ/ixvqVSiV+9dufV10zqJfWvMSP/9+/0dHRQWdnZ2j5RYS+vj4eefxWurq66OrqCoWiDzuAUqlENputOu/HccFXaTVGtT0lBaHbdhgMqx7mZCpqsLp7/fBrtjsAxlKmZd+uEpFP4K6h+AJwJe4q27GrZrf3XU1iBR/44IDnEfBUhgU7K0wgBSFjuslBSzHqoS9vVFurYTD0yxQEX5SV6M9hjmp3Ct4DVAOzEpzxw5ySzFKJ+5OyqlTUdxLlCY2DelX1NOAJVb2DmEUdfFkLcYwV5kaHW0wuIPzFEkzQZDKZKoik0+mye+W7p0nKYn4Gz5kuPISvgRcsiwnDMFCqUgH4KMvQl7lwQy1FLdAwcj/h1mSS+ok6nsqaAHXhNwT7P4JEQwcsEFugqHaoMAstCAQXLJXWlWmFVVtc0cNXwv6pzeOZM2cyNDRUcS4szdy5cxkaql4TMQqCUWG+hegv8BAHQzesutzRg9yrpzbW+mNo1E0P5hFMMwFgMWq1u8sMlETkD0CPiByL+8KpWLX9XU0UJbUmwtxQf38EJpWusmklHn74UeXZIKbFaLrTn/j7/+Q9b/t0TXc5m83yT2/+Ax9/+y+YOXNm5H0tWrSYbN+/cPSBf18FQxN4MucI/rrjBKbPmB0JyGw2y/auGeQ7p9V0l4NhSVww87S52IUZNmIsRH9fYZ9TAXD1KIm7PF51JiI/FpF3qerngW8CvwW+r6r/niS9BWILZP4WoiwQf+vq6iaT3YV0urIzwodgJpMhk8mwcOEi3vzqryOrvsYbzvxIqPWVSqV47ZlvpvOxD7Jf7kvss/c+kT/Md7z5n/jLT+eQKsxjYGCg8qRhgO27y5t59D4o6fZQYPmwW7ZuFnetyJBOUS57EOhdXV2ct7KL6wdmRE4DFBEGU91ctnIuG5mBOewm+KkIz6yZjaQ6ay56oaqoI4jU5xBNNauvHrUrEIHNqnoVgKreoqrfUNVbkya2QGyholzkkdkoKX505T588cuL6B/YKdTV7Onp4YQj3sG0jf/OLT/Zk533LOGwowxKEya777Y7Z8z7CTs2w/pFl7HqWXdee9C9WbBgAbM3vZfp8+BPz3+OQsEd4lKGiduXwaxZs1j51/3YY394ceOtoUBMp9N0d3ez7MU0h+0u5HLDVTD04wF0pISchrdH+tv0jhS/fU65f1OHV2/hD1iuOJOvXT2dx56dVy5/0DpUVcTp5erPHMjmZ5dUNWfU+0Bbi9FVGwMxtLFZRC5KktgCcQyk6rcrBn80cMpJQ3z4H/qZPXtHiCuaYc9dD+MPVxxBbghe83eP8ljfe7nt7uvKVlc2665S3dPTwyfechNP/rGD/f7mQb71o8pXkJo/0Pe/5Vs8+acMe5/1KPc9eHegrCO/pyMOfBNPLhUOOeEFtm7tq3DvzW3ajHmseFk5YFEu1DI0LcGOFBSMjpXQQdylPMcuSvPiQMrr2KF8D+Z9dGe38ZZXKz2dw5EwVFVS2SHO+cBmZu26mji3OWjRR7VxTnUo1vpD87dx0vtE5NnA9hxwWZLEtlNlFHKtvVrLfo30kqZSUrZOUinBcdwHbY/dV5PPF8jlho2hLCMdK8+9+BjnfGgJK1b/gdvufZKBgQEy6UwFbNLpNI7jsGnwKY571zQ+/f3zymHm0BuAOXPm0NV/LHsf43D17z8W0s7mfqbTaUo7jmHBTrB+x41V//7mD3+IOUzvhr3m7iDtRA+9ERE+tCTNnl1a8cCUXWBjgdf/PLofLfrTGavL529nHPE8qk4kDFWVQqHAzJ1fRlJSkZ9/7TCLZqoDr5ba3Er+JdXwExIu7mCB2ASJCI5jLjIquEva+/DTikUbqqEiFRBx3eE02awLvkdW/JLt27e7izJks2Q7smSzWTo6OsqvIR0eHuY/rrqA7u5uduzYUQFC85pbtmzhsv87mmOPOpkXX3qRjo4Or9yVsCmVSixf8zmOOPtkHl/xIj09PVX5+eUvbXuar79uPWlNI6lUFQhNi+HN0wY8qy98YHT5D8YpglT32IctruvWfTgMR/IeySPYAx/l7lmLMFptXB+bVPXxYKC4L66PlXWZmyDzofVdvLAHztyPcj9dV3OkMyWTyZThl81myy6yv1+2xFJumoGBgdhe2s1bNnP7Xb8zyl89XEVV2bJlMw88fHOoFWzmm06nKBTyOE4pkSsV1qNcb1tUcHXxMOux+p7c7yfqPsJgGNeeOFXVxi7zcSJStZCDqg4nSWwtxCZo5Lmo7ULHQdF0M00Ympv/YijzuPzCKbS8Eo6/EnZwCqHpngYhYr7RLxgWvSJ3tJUVBkGzHTVOUauKV8Ouuuxm+pF6r/wuguUIg2HY91dVA1MQjO16z6p68mjSWyA2oLi2w7C4Iw+f2ctcPQUvDIpBMPpQLBQKZLNZCoUCmUyGUqlEJuN+pb7ba5YjTEG4+G/yM8ES94qCKKsqbjPLUKtO3bbYcCsw+ArXYDyz/v28apXH/J5q3dtU1mSuA+syJ9RofgDBh9GfEhc10DkOhP7W2dlJR9b79MLMtOawHP86pqLAUiqVKl53ap43w0yX1c+vVn3V8yC5+QOB2SfV762uhmEQipXXrwfYI39eVpWq5w8vJp9jRGSlxAyLEZELRORhEXlIRC6TFn0pFogRGm19h7ljbudJtVuZBIomCMttiplqFzqbzdLd3c2bz72gnE/Y3GGAY458JSccdWYkUIrFIo7jsN9eRwCUIelDMwpCED29LqqTJMod9l/FGn7OCXyGtx2acv8U4h9k/6ur93cwVeDZjDZEcd+b/I/Atph4B+P2HL8GOBY4EviHmDQzRaRbRDpF5FwRmZXkvqzLTH0ucFRaM4/gPoDjKI4xHtF9KFMU0jMoogj9pEul0DbEWbNmsWjRzhSLixkYWEgms4k77vguxWKRYqlIyRmx6A456AjO3O8yHrtpLue9todf3vSDUMvtjFPO4eDBb1MaVNYteT3r1q+rssBUlaMPeTXP/uItHPGG9Ty1/Yoq69Hv5V696nBEShx4yLMVefiQGSmDg9/WaoYHQTcwlEFEyaaHvbJUWoW+az+8o5ueGQVy+aFyb3PwXkUEdSCdzSBC5J9SnBsd/IObymrS/T+oqr8SkTtj4l0M3KSqm7xr/xj4JPDdGmkuB76CO9xmIe4Lp94TVyBrISZUrTY0U1E/FBEhnRqxEt25yinefPU03vrTXh5cP7vCQkyl0mQyWebNO5H7738LV191NldfdTi/uGEx27btQ2dnZ8U2a9YsPnD+N9h320944Pq5HHpWH+v7nim7zeYD/ZZz3sM+G7/D+pUpeo+9ledfeL4KhqVSicMOOoGXb3kzIpDrfaAMQx+CvgWZyWR48P4s69d1VoTHubdR7ng63cFHLpvNfY/PMsBbDUNxpvGTS/bkuUd2q7ISTaUkw72XH8jwppmGKxwPQ9+9bvC/ctIqeZNDbanqSwkveQywwjh+AjhIRGq9XewB4HHgeFX9W+DpJBeyFmJCjd6FhqCrpgpff2MJUWVR9yDp4oh1mM1mvGEsA+yyxOHII/vpnbaBYmkla9c+Wx6DWCwWKRaKHHrQK3jiulNZsIey7/k38J/Xf4XBwcEyEH1g7L777uy++d/Y0ifMfcPVfPMn/0Ymk6loOyyVSogIC3PvYdlm4ciLbueh5X+kp6eHfD4f0qGTJZ1Rluyao1gsIuK+BsCXqpYXnjA7OMzz/qdbjiE+9vYOFs4e8CBYCrUSYTt/80999Mx7IWAdVl5HUsrxF79IR+8wkDb+lKLXjvQ/g197sOxT1VJMOKxmnog8ZBxfoapXNHC5hVS61Vtx3ax5wIsRafYF/hW4xzuen+RCFogNKs7NjnroR865g7cPmbN1BGxGO6Lfazw8/Cy77f44g4ODbNw0xMDAAP39/XR1deE4DoVCgVJXiUeW38trz/8df156PU/91v0z9Xuci8ViGYovvvgi6w7+Lzp2g8t/9u2KXumyy1pyyOfy3Lfu0xz79rP543030NPTU7YMC4UC+Xy+PCZy06aNvOo1t9PT00OhkKmyEPy1Gv0yBM8Hgeg4Dnst3ggopZJpXVZbnD3zny23TbrXCMDQ+54yPYOI+K8qGFmCbKT5InpQdnDs4liqXYGbsFybVPXoJl0yesxXuK4B3gp8RUReDzyX5CIWiDUUhF4UBMPaEcPimHHdh8x13cIWTAi2I3Z0dJRh1N3dXeG+qiqOOvzurh8wNDREV1dX+Vrl2SRefFXl+pu/BVT+y5eB6LhAKjklXl63mj9s/Qnd3d0Ui0UKhQK5XK5ctnw+XzErRUTKryQwXeFgHP+6fp2YQKx0fU032QVisEe8sje5cohTsJ3QrdsU4kGwfBzTdmiWNbg/VTXGdbABmGUcz8IF5MYaaR4HHgFKuJz7aZILWSC2QLVAGtXeYsJQVclkMuX2uWw2W9FuVyqV6OzsLHdqBNvnwBuDFxg07VtoQBlUYZaso06F+1wsFssWoW/N5vP5QA/2CFiCQ2R8a9d8/YFfLvO6wd7lEddYy1A0B2GbMHTzcqEYHOsZrGd3vcnwxXVrtYHFQWCqgDJpG2ET9SCwn3F8IPC4qg5FxIcGO1UsED3FucCN5gmVb04LWjCpVKr8cJtAHGlLzFaAzweib+0FB1KXr6W470VJCfl8vsJK9Mvhl7H8A9dKKy3YiWJahOWFYY32UH/MoA+ybDZbzsMcfF7rgTKh6FqFIzNRSiVFpHpYjWmh++1+YdMDg2GVEI/uWY4qaz1QmGywbOXUPBGZB/wKeKOqbgF+CNwqInOBPuAiXODVktmpcqyIVL+bN0STHoitAF2tvJO2LdayEn14+NZVsVgqT9Ezp+qZn6FWU8oFop93sVgs5+uDsQJQXpOaCTUfiIVCoWoWjVlmtx2uMq2ZPuw9MGGQCFqJ5sDs4CIOvswOEb9uzU6TcCiOxAm+O6YWvCYb2BpVM+pBRI7CHV94OPBpEXmDqp4HdAP7Az3AFlVdLiL/DNyCO27rbuB/YrK3nSpRaiYU49oRw46DrqEfFmYl+hbiiJXmkM1WdiiEDV0JdiTgLcIq5bay9EjnjQfGUCvRUBjYzLZE0/31Oz38Mvnx/QUoggvFmu56WP34HSVhAKz13Yy4wJVTImstMjFy7+Hth8EyTnU1y2VW1YeBU0PCXyQAMFX9GfCzOrK3nSq11CgURwvTMFBCtetsgrHSQkyjmikPPSkWi+UhN5GWoXeNMhDTla8wNa1EE4omJILtgD5ETWuxctEGytc3IerPr/bTR8OoWn7PsXnabRsMB1bwHnwoxoHQtxjNfILfV9j3OpU1jqvZJJKq/kVEngYWAXeqan+SdFMGiEnUCPxquc1xMPTDTEvLtBL9/dXOHvT05piWGmRmaQfbt26q6DzJZrPsttvuOE4HpVIPTqmHkjMNdbp44skfsWrVKtLpNIVCgVQqxezZszn0gFewa/fprNm6nBv/dHlF+XxwAbzioDMhXWLZyntxSiNg9FfT8dP0ZOYxu3MRW4ZXl2HotyH6FmJPTw/rH3sFc5cMIjMfqbjvKBD5ABypVyrhVcqSyaZBnJqW4Jbts1g8f7jcvlirIyUIO+tKV6vd71dE3gF8BlgG/ExEDlbVr8Slm1JAbGV7Yq38TRCax/5+9Tbi7pHK8OlfdwAdoNOBhXzhpLn0ZP9KqeS2q3V3d/PH29/Eyy9Xl+n0V+3O2rVryefypFNpUpLijft/m79esxNbHDjkLUvIZn8UWZ4Fqy5EegZ5TP8yYi0WiuW8BGH27DkUfvNeOLRIYe9vByzLEtms28vc09PDI3f0cvipWaZ35SLB1NHRQf/S/Zi++w7yvS+U2/jMuvM/tz22O3OOfp5UutIidNs6Xdd5OD+NS74xnQ+8s5NjDt0RYknWHnJjVa0JUDeHqeoBIvIpVf2NiBySJNGUAmKzNBqwhsHRtBIdR0mlQDXlgaXEr9+2hS3DGbbkMmwZzrBzxw6GujpIpcA1sJSzzupjOJcnJf3ADgrFreRym1i5ci09PT3uaweybq/11o4HOfb9h7ApdS/3rbyDjo6OKkvWB8bgQf/HUG4QVgHitSuqO07RH/C9adNGDn3LnawdfJz8cI6OjpFe6my2w3Ob3VW99zlykFmL+xnyer6DYPJ7wzMHPMqwCDgu6MyymXHnH7EJMv4rW32X3N33w2Z0Ffncx3Lssrg/9HpBCMa5zBMABi3VBPmz8Ge2+D/sdFREUxaICRQHwKRuc3DfjOuPEUynU7jfofv9qRaRQj8ztci0VJGFmTz54XzF0v+5XI5t239PPp9jeHiY4eEcQ0NDFAp5CoU8PT09rkWXSpFJZ7h3xf+Sz19ddm07OjvdeTOBoSwiwt2P/aqiM8SX7wr7YFy29mY6OrJkMtmyizvSrpkBlC1bttC55y/YUkqRLVb2OpttlyOgSiHiLoyRTldbcalUCs0MlO/NHV+YLo8zHBl4DbvttN27h/BxoFFuctRxPZoA8Khb7d6GCCwUkcuBxSLyDcJnulRpygFxPN3mqF7nuPZEU2F5q44sjy8iZXcxlxt5t3M6nSbbkSWfz1fMRy6VSmS92S7mwG6zPCISuoyYXx53c7w50CmKxWrAjIy19Oc1UzHu0u9MAsp/DiPLpVW78mEL6voWYXWYC0p/yl6wMygKjElB1irgtTNI27lsnj6Gu0rOocBTwI9qxvY05YAYp6TATBIvzkr0j819H4g+IILDcczeZPfTTzvS2RAcf2dOtfMHe5urbJtLepnXCJYpOIwlCEZ/dkmp5JBO+4A0ISP4fPeDTCgF63NkfGD4mpG1tuqOFdM9jht07c90Cf+OrNq/HlS1BJQXkhCR04A74tJNSSA2YiXW4zbHWYlRPybTdTatxDCNgHYEHP7D789M8eEQfDGVvziDPzfaXyQiOK7RLJfZ3hYET5iF5U8NNO9BvVknUc+SCUdzmp8J4iQwNOOGQZGQBR3Msrsf7f3Aj6dML6ZdJSJnAx8AenG/zF2BveLSTUkgNkuNWIlx7Y3+cRCIwU8/HkA+X3kcBJc/5c63Cv2lu8yB2hUvq9LwwdAm/ILgGenVNdv/qoerRHVSBC3Ramu0UcswfPxhcJHYYPnCrEOrEU2AevkMrtu8EReI70qSyAIxRI22MyaxEs24voJWo28lhv0TB9sU/bzMB94fb+i3+/mr5phrGBYKhQogmq7zSLtguKUYbMczoegubFvZljcCt3QV8JJbe+67qpNYhGFgDLMGfQvbbKMMfjdR31mtsKmgCXDfS1W1vBajiFyTJNGUBWIr3Oa4dFFtilFpkrolIkImk0XVqXovswmITCZDZ2dneeqdCUN/5kvNudGBTgjzGtOmTSOXyxlzlivnPafTaXp7pyN6CtNnbGf7jr/GgtBPZ4I2ujMlepqe6e6HW4SVbYaV50anCQCOhjQB7muRiPyMkZWyTwLOiEs0ZYHYLIV1ktSyEqPOB39gpvsbds3g9lJqPl99ahqvW6KcMGcz3dvXMjAwUPGyqmKxSM/0A/nriv3Yc5dtzJ72HH2bn2LLli1VnSv+4g+7pd5Gz64v8/izd1eVzQfNznMPZe1dZ3HouQ+w8qX7vPY/vzNnBFrTpk3jxl8t5LgTZjFj9rJI+KXTaUr9h6JdG8hM31YFxCQWYS0Y+hCEyil7Yd9D8J6tJkYbIu6UvR8ax3ZgdpxqWXxxHSOjvVaw3TDMNTWP436A3Wlh75nCd55w+HppDkfPn8d79sgxd8NjZYgUi0WynbPY0ifc+9AsHD2CJYuP4A1nrOfl1f9XsVy/4zgsmLuElT/cl2z3Puz0lmfp27q5qowiQs8MZdv6NE/f+AoWnLGKgYH+MhRd99l0qQGV8qIPYVt3xy5894szOfGsXo48Y1mkJRh02YO961EwNN1l/17Mz+B+re+0VWp3ALd7+YAPqOpyAHGXE/tjkkRTGojNUlIrManrHNWmWMtaXFTq4xM7beHju3Xz4MBM/rAuzUARdurspOi9q6VUKlHMLeO1Jz7Bm8+cz8ZtO/Hsi7PJyHZ6e6dRKjnG8BuHHYNbOO6jS9myJkVfYZju7u6KMvnX37xtDae9/yVW3jOXbKqHrq6i4TaPwCqbzZJKgao7PS/M7c1kMmhqI3/zwRnstn8fJSdbl0VY2XFS+216SV8NMAEe/jHXBKiT84Dl3n4W+DIJOlYsEGuoUWuw0evU+pH5YxLNNMGOAMdxoDDEUakdHLWTO0ZxuKODTCZNPp8ik3HKYBzoX0+HrmHfnRwKw0pXV3cZhP5y/aWSw9odj6HTlS7tRqR6PrZfhk3FPzH9aIcS0NnZhTkA2ocYwGvP2ci03iEc7ahqD/Q7XTKZNDvvswZHUxXLhiVtJ4wDYbD85mdwfzSaANBoSO3sMovIobjrKx4uIu/0glPAjCTppzwQm9W5MlorMQkUYcQyC7u+uZWn5XVkcRwXNv5SXGZbYbFYMuL5FmJw+X5/vcXKl80Hl8zq6BgZ52eO/TPd2WnTn3CtxXQnIpVjGoMudtjwm0ZBGAS4v29+BvfDjuPCp4ra+P5nA3sYn+C+V+WbSRJPeSCOleqBoh8nSqa1GPXw+y626W47Trqi06RYLJHJjFiEwTfajbx/Jbhkf+W6hKq+1RCEUnCqXOVxLeglcYuTgrAWDNtR7Vw2X+1aRlW9C7hLRK5V1Wf8cBGxizskVS0rMepcvVZi3LWDFmKYtejH8UEVBUYfgiMgHHGV0+m08fKn8Nd7BmesmDB0HN9SxCvDyLQ8H4omEM351bXglhSAURAMG1YTdmx+1toPO44Ln0qaAHWwUkSOAKZ7xxfivnCqpiwQx1D1uNEQDkUYAaEPLv/Yf/jNfROQ/qc/+Dqdrn6tp/lpLkJrfobdl7cHgPnO48rZK7Vf+pQEgEnbCc19s4xJYThaTQBgNCy//ttcv/U+N3mfdthNs9RMKzEKisHzcVD0ZcLRj+vPcvF/uEHopdMj723xxxuGWYR+Xv45X2HTB33L0x/bVwmpkRc7+W8RrBd8SdziOPe4HgBOZqA1QxOgfrap6jv8A89ajJUFoqck7m2jeSWBYlSYqSg4+vHDZpaYcPPheP/sHg7JOWS37SjDMWgVSjrDtZ0LOLHDYafNa8p5hunx3iU8kOvg7/SFChirVlqLqVSKH8+czcJ0ijdu31oTenEucRwA64FhUlc57txU0gSoh4dEZE9VfdY7PgxYGpfIAjGh6rES60lvnkviRoelC+YbhJyZj9PZwc/zO7gO5YtzZjFje39FPPBgmu3giWHlui1Frl20iJk7tkSWf1AyfOfFPOccsIglQ5vKZSOwokwmk+HpfJ6uDnfVnTjgRUHQzz8pBOuBYTM0mjwnAGiA9l0gVkS2AFtx228+KiLq7c8AroxL3553NU5q1I2KgtVo86r1gIdtYQOTq+YHF4r8R3om01MpPlnsY8Ps6eWpc/4yYZlMhpRT4tIZeRZ1pPmHTUJp+qzy8mHB7VXpbew/Lc1lG7J0dnYa5zJksyNTBzs6OthWKjEjZOWaWgOs6wVl2GdYndf7vU0UWLVaUd9L2PcxDvqQqu6pqnt4256qugfw/iSJLRBbqHpcs7j9OMsoCMKoDot0Ok1vvshXZAZ7ZDr5bG4LMr23amWcTCZDVzHP5QuEIUf5xLYs3T3TKs77EE2h/MsuaW7ZXOCp7OyKc26+I8dbSw4zUlIFxKgVbJIAMg6KSb6LqLB6NVWg2a5AVNVrvfK9PhB+Q5L01mUOKM4FbrbrHOYm+/tAqNtcq6MlLD//03SJO4slPp/pZWX3dLIDOUinK877n3OKw/xgSS8bS6DFITKZ8J/M6ZlBLtljGktSOTKSKV/fvfaIC/2qnm72yaRJ40RavMEw89ism6QAbIabPBYP+ESC6QQo67dF5J+MY8V9VdoXVXVNVCILxBAlhVuSdEnCokBonqsFyFplCX6aeacch32HHMSAoXne39+rNOQuNZzJlAdkV0v5u+n93rjEdCisHKfEJ7NpUk4JDIvPPx8Hv1rgSwrDWvXVyLmpqnZtQzR0HXA78CywJ3As8GfgEmq4z5MeiI3CrZl51gvF4HGU1QfJARlMa45V9BUFxbDjYL4j8SoXTTCBlUrFW3pJrcAkAEzqFjcLeFMFnOPpEtehIVW93dt/TkROVNU/icgraiWa9ECEaqsraZpmWYlR4fVCEQiFYT1lCgOqr/KQm1H+iSSBV6OfcWFRx1FhSe+l1ZoAgKnQBCjvMSJyNLAS2Bc4VtwXfNd8r8qUAKKvZlqLtfJqFRSDx6OBZFh+acN1bvQHHwWmJFBMeq7WflTZo+7HWoeNqRkus4h0AZcD++Oy6F9V9ZaQeBfhvh9lqxH8WVX9c43sPw98HzgIeBz3hVOHAvfUKtOUAiLUB8W4uGMJRaheeqtRMCaBrKla91hLSaHW6Pkkx7XKmbT8VtVqUt18ARBVPU5E9gXuE5EDVHV9SNyPqeqdSTNW1aW47YYAiMhiVV0LPFIr3ZQDIkwcSxGqp8klBaX5g62nfTHqXFIFy1ALUPVaf2MJw3o02rwmGnijvtc680gB78VdyBVVfVpElgLvAC4bRb6HqeqjMrIWoq9zgLfGpW/7rqLxVqt+rEkf3rCHPi7MPzbBExcmUt984bC0cenD4hDyfuTgPUXdXyP1WSu83jhTWUl+EzHaE5gLrDDCngCOjoj/fhG5S0TuFJF/qJHvR7zPd+Ouh+hvc2JviilqIcL4u85R56LCIHxRhaRhZnjYj7UVbYdhlmmrj+PCk6gRi3iqKWEb4jwRecg4vkJVr/D2F3qf24zzW4EDQ/JZD9yCO/VuPnCniKRV9dvBiKp6sbf7OWCZqm4DEJGDkhR4ygIRxt91jjoXBrao8CTAi4NdvR0xUYqCVbDcSS3hWsdRYbXC484lOd9sTVSgJiz3JlWNsvh8BR+MqoxV9ffG4QYR+TauJVgFRENXA6/CA66qPh5fXOsyJ9Zo3axGzjUSHgWOWq5Mo25yVL5R4b67bHJ+LN3hZsNnosJstKr39xChDd7nLCNslhFeSy8Au8XEuVFVnzPKfFqCfKe2hQjNdZ3j4jRyLs5ajDuX5HxQSe4xiWqDOz5us/4o4s7VE6eZmshAbcKwm1XAFmA/RiB4IHBTMKKIfFJVv2YELQRejsl/kYj8L/Ckd3wScEdcoayFSPPbjEZjKdZr/cWdM88n/Odu+J8/Kl6wx9nMopbVWs99jjUMJzLMmqHRWoiq6gA/AN7j5bcP7tvyfiYiB4jI7TLyHpSzReQUL14P7tS7a2KKuBj4PbDa27Ymua8pbyGOh6KsPvN83BCYWtakr2YMoxmNakOxec0EScswmjijid+qPMZTTSr/F4DLReQ+XBa9XVXXicjuuIO1s7hvzLsMuETcoTq9uHOUvxqT97tVdaVR3t/XiuzLAtHTWLrOSeIkgWat82YcX83qQEp6vWD4aCznetMkOW/VmESa804VVR0GLgoJvw/Y2Ti+EbgxYdkywFnARhFZDXwD11r8PAnaJy0QDbUCihAPrSRQS2Ltjbbtr1lth/WkaTYkk56vN16j8VuVx3irje/hWmAe7hjHR3DbEJcCXwLeFpfYAjGgeqEYtT6gVbWGh4fL+62EXTvDcLKojetis6q+zXOvr1XVSwFEJK5XGrBAHLWKxWJigCaJ18y8mpEmTu3SIVVvecbrgW5jkNSlNr6PteB22ojIX43wQpLEFoghqsdKrCd+s9zseuKFpRlrNQNk4w3DNgbAmKtZbYgt0mtEpNfbP1ZE5nn7xxHfEWOBGKXxhKIfD+pv12t1x0lSRQEkOFulWaCbCDCcTFBt43vJAwPevjnu0FqIo1UroQjJ4FWvJThegKwHSEmnCjYbctYybJ7auF4+qaoPBgNF5KgkiS0QY9QqKDYSF+oHXJylVq9G8yCMl1VYb9zRpBmLvMZb7ewyh8HQC384SXoLxARqNRShfgtwtJbfWD6grQBXq2HYTI339VuhyXhPYIGYWK2Eoh8f2t81TqpWDmkZCxhO1ge+WZqs9TPpgdgsi8rPq5VQbDSNn87XeMFxLOAz0WA4WcExWe9r0gPRVzNdzVZafmaaetOFpTfVzLUfxzqPsRpEbWEYr3ZuQxytpgwQfTVqgY02j9Fafs2ycMdbY9G50Q4wnOyarHU1OTEfI5GmvCSn6df90Y9+xMknn8wpp5zCEUccwd133w24U94uvvhiTjzxRI477jhuvfXWhss9HvLvuxGwTVQYTlZg+DK/06htImrKWYimRmstjsbqC6a74YYbuOOOO7jjjjtIp9NceeWVrFu3DoAvfelLqCr33HMPTz/9NCeeeCLLli1jwYIFDZd9LNSMITpjcT0Lw/o1WV1mabfeyWZLE9zgaOtgNOn9tIcffjjXXXcdBxxwQMV5x3HYaaeduP766znppJMAOPPMMzn77LP5x3/8x6aVoxkaj7bF0V63FfCaQEBsqKCHH3643nJL1fvkq7Rw4cKHE7xTpa00OTFfp8bDfTbTbty4kRUrVvDoo49yxhlncOqpp3LFFe7LyZ599lk2b97MfvvtV05zwAEH8PDD1eNMfVfl0ksv5YMf/GA5vK+vj56eHoaHh/nyl7/MmWeeyZlnnsm5557Lyy+7K7E7jsOHPvQhTj31VE4//XTe//73MzAwUHWNsOs1qwnCwnDiyLrMk1zNcJ+hMStt9erVqCq/+c1vuPnmm9mwYQMnnHACM2fOZMmSJQDMnDmzHH/WrFk8+eSTUdmxfPlyTj755HKZHnvsMfbbbz+6u7uZPXs2t9xyCyLCVVddxb/+679y1VVXccstt/D8889z1113AfCWt7yFTZs20dvbG3mdZmi8gGZhODpN1nu1FqKh8XL5crkcjuPwwQ9+kEwmw0477cQFF1zAlVdeGZlvLfAuX76cww47rHz8yCOPcMghhwCwyy67cMYZZ3DaaafxrW99i6VLlwIuZJcvX85tt92G4zj89Kc/Zdddd637XpJqNFbEaC0QC8PRK5VKxW4TUROz1C3UeEBx9uzZACxcuLActmTJEtasWVMO27p1a/nc1q1bIztU8vk8q1at4tBDDy2HPfbYYxx++OE888wzvP3tb+fSSy/ljjvu4Jvf/CaDg4MAHH/88Xzve9/j61//OnvttRff/OY3m94m2Qx3ajybN8Yyz3ZWEnd5otaJBWKImgXFpPnss88+9PT0sGHDyCsfNm7cyOLFi9lrr72YM2cOTz/9dDm/J598kqOPDm+rfvLJJ9l5553p6ekBXEvy7rvv5tBDD+WRRx5hxowZHHPMMQAUCiMrIm3bto1TTjmFm2++mdtuu41rrrmGa66Je7FZMo1nG2MwD6vmyAJxiqlZX2iSH0dnZycXXnghP/nJTwAYGBjg+uuv553vfCepVIqLL7647D6vXLmSRx99lPPPPz80r2XLlrFhwwZWrVrF0NAQn//853n++efZbbfd2Guvvejr6+Ppp58G4Oabby6n+/Wvf80PfvADAPbaay923nlnSqXSqO65WQ9GO8C0Vt5TUdZlnoJq5o89Lq+vfe1rFItFjj32WM444wze+c53csEFFwDw+c9/HlXlhBNO4B3veAc//elPWbRoUSh4li9fzplnnskZZ5zB/vvvz/Tp01myZAlf/epXOfLII/n0pz/Na1/7Ws477zz6+/tZt24dF110Eccddxy33347r371qzn++OPZc889ufDCCxPfWyssg3a3CqcqDCezy2zHISbLoxlFaVl+pl73utdx8cUXc95557XkemPxQ2+mdd4qTdQHPqCGbuLII49UfxZVLU2fPn3CjUO0w24SaLRDcsLyg9aAcfny5ey///6h12t3jaVF3q55TxRN1jqwQEyoVkCs2Xn29fWxYcMG9tlnn6bkN1Zq9sNlYdh6TdQ2wjhZILaBmgXG2bNnMzQ01IwijYkm2hAYC0NXE7mNME4WiHWqle5uK/NuJ000EFpVa7LWtwVig2p2u+JY5T3emqgwnKwAaFTWZbaqUquhCJPHWpzI4wAtDKs1WevEArHNNRnAOFkfnqkq24ZoFamxcm8nohvd6odmsj6UE0GTte4tEK2samiyPvij1WRtQ5z0M1WsrKyaKxH5AzAvQdRNqnpWq8vTTFkgWllZWXmanHavlZWVVQOyQLSysrLyZIFoZWVl5ckC0crKysqTBaKVlZWVJwtEKysrK08WiFZWVlaemjJTZW85SwdlEwqoj1gB9baqYy+s1rGKup9GWv/Y3B85p5V5mNfEv4ZWnPfTlc9LORniXV+8AME8r4gZXxTjkkZ8I9xIG3qe8DCv9IiMHI+cVyNe3HFIOo2IoyFpVCs/ATTkfEi4nyeo/zVV5Fe+ZiAP/CGy/ljZYJiXp/uj8OOE7FfEC+wH8w67RvA6kfG8AMeIDzwMN0+0AcpTVU0B4iCbuDj7EE4aip3gpEe2Ygc4acXJmMf+vgbiasQ5LedTSkMp654rpUHTSjELpbSiaXBSSqFD0ZQbv5AFTXl5pZRiVt3ypHD3U+5+Ka1IRkmllHQK0lmHdApSKUVSSiajpEURLyybcUilXBhmMw4pAUm5cTJpx40nSjbtIKKkU0pKIJ1ySIkbL512SIuSQt0wL31KFEHJphxSqHsNGdlPo2SkhOCmzYqD4IanUDI43jk3TVbdzxRKWkuk1IurSkZL7jk1zql3znHPiSrZknsupY53zkEcxwvzzjvu+XTJIeU4bj6OQ8Y7TjlKuujGS6t7nPGOxXHIFrw8vLgUSuAoOA6UFIolKDluWMHYLzkjcf39krdfMuJW7HvHeeNcyfGOzf0a54rBc95xrjKdJJvVYdUGsi6zlZWVlScLRCsrKytPFohWVlZWniwQraysrDxZIFpZWVl5skC0srKy8mSBaGVlZeXJAtHKysrKU1NWzBaR5cDw6IvTMs0DNo13IWrIlm90avfydanqweNdCKt4NeslU8OqenST8mq6ROQhW77GZcs3OonIQ+NdBqtksi6zlZWVlScLRCsrKytPzQLiFU3Kp1Wy5RudbPlGp3Yvn5Un+xpSKysrK0/WZbaysrLyZIFoZWVl5SkREEWkS0SuFJH7ROQhETmzRtyPi8jD3vaJ5hV19OUTkXNE5BYR+aOI/FVEPtJO5TPiZ0XkGRH5QruVT0ROEJHbReRuEXlCRD7cLuUTkekico2IPCgiD3hpeltdPu/ax4jIShG5KCbeBd6z8ZCIXCbir6Vu1RZS1dgNuBS4ytvfF9gCLAyJdxawAujythXA65JcYzRbHeVbCRzh7S8CNgNvbJfyGfE/CGwDvtDqstVZf3sADwDzvOODgMvaqHxfA+4G0rh/9ncAl45B+d4EXAs8DFxUI97BwDrcgeR++T44Ft+x3ZJtsRaiiKSA9wI/8gD6NLAUeEdI9PcD16rqsKoOAz8D/j7uGqNRneX7rqou9eKtw/1B1rTWxrh8eBbN24DftrJcDZbvn4CfqOomL+7jqvrxNirfQcADqlpSVQcX3ke0snyeHlTV84EdMfEuBm5S1U1e+X5Mi58Pq/qUxGXeE5iLa+35egIImxlwTMJ4zVTi8qnqfwaCuoCNrSsaUF/9AXwc+DZQanG5fNVTvlcBHSJyk4jcIyLfEJGuNirfTcDpIjJNRHqA04H7W1w+VPWlhFHDno+DRKS7+aWyakRJgLjQ+9xmhG0FFkTETRKvmaqnfGWJyAzcH+hPWlOsshKXT0TmAyer6g0tLpOpeupvd1yL5p3AqcCBwH+1rGSuEpdPVb8L3Ak8C6wG/gp8qaWlq09hz4dgX0LVNqqnlzk4YDGqMXi8BjYmLZ+v/wC+pKrPt6g8QSUp3+eAr4xBWcKUpHyduE0im1S1gAvDizy3ttWKLZ+IfAY4DNgN2BXXuryo5SWrT2HPh+1YaRMl+SFv8D5nGWGzjPBg3GC8Vruk9ZQPABF5H1DwLIpWK1H5RGRPYA9V/eMYlMlUPfXXB6w3jl/ChWQrLZx6yvcR4PtGG/blwL+1sGz1Kuz5UFr/jFglVBIgrsLt1dvPCDsQeDAk7oMJ4zVT9ZQPEXkTcBrwUe94nzYp3ynATiJyp4jcidtjf5F3vEcblA/gESpd1flAHre3vlWqp3wdQME4LgDTW1e0uhX2fDyuqkPjVB6roJJ0ReMOe/iJt78P7gOwCDgAuB1Ie+fOwm0o9ofdPMnYDbtJUr6TcXuW5wK93nZVu5QvkOZKxnbYTZL6+xvgcaDHO74K+FEble964H9xXVDx9v9vLOrQu/6dGMNucC3nPwFzvOODgbXe7y/lld0Ou2mjLekX3eU9oPcBDwFneuHHAWtwF8D0434cdzzWw8AnxuQmEpbP29fAdme7lM8L6/AerHW4HQO/bLPy/QvwGHAP7rCR6e1SPlzr9f/h9iw/ANwALB6D8h3lfWdbcXuRf+mF74LrDi8x4l7gPRsPApfhrSdgt/bY7OIOVlZWVp7sXGYrKysrTxaIVlZWVp4sEK2srKw8WSBaWVlZebJAtLKysvJkgWhlZWXlyQLRysrKypMFopWVlZWnhoEoIq8WkUdEREXkLn8OrjEXt+Xy1uU7tcXXuERE1kUt5x9SD3PGo5xJFFLWP4vI40lepRBXD3WUoSV10co6FpEvjOa+R5veauyUaTShqt4qIh/DnRv8KlUt+udaAUQvzytV9Uoj+G+JX6V4VFLVL3kr0USdj6wHQy0vZxKFlVVEDgKWisiTqnprjbQ166EONaUuQn4PbVHHVhNbrXKZP9mifCukqtt1Asw9bOdyqurjwDLchTnG4notqYt2rmOriaOmAlFEdheRK1X1ARF5n4isFpH/FZEfeK7anV68S7w33/1RRG4UkZ2MPPYSkd97Lt29IvJFEfkqcDjwac8lf52IfCLowonIPiLyB++NcH8RkbO9cLMs3xf3jXs3mcvf1yrTKOukopxxZfHu/xbv/v8kIickKWOt+k6gLN6yWVF1GHFvoeWJKktIXVxknLtTRDaLyE0J7jX4e/h8M38LcRKRi0XkZhH5noj0icjTInKgiHxURF4QkU0icl7S/KzaSKNZGQJ3GXkF7sJd7eM+XDfGP/8F3FVb5uPC9z+88A9DeWGJi4BrvP007vJh7/KOZwIveft3EnijGcYSWbju/wo/DrA3sB3YyyjLGmC2V5blwNuNvELLFLxOTD1kIs5XpI8qi3cPTwLv8eIdCmzCW1GmVhlr1XetsnrHReAVCeoweB+16izquze/s/OBE739g3GX1z8k4b1W/B6a+VsIqbMvBO77W7hLkL0K9zf7C+B53DU2s97nX6PS2619t4bbEAPy26N29758U/eqqr8i8Ke8zxeBO8Rden4G7pJX4C7ntDfu2/pQ1W0i8jcJy/AK3CXjf+qlXSki9+Mut+S/V+N+Ve0DEJHluK/V9BVVplYprCyvAPYCrvHu4TERWQO8HrguYRnD6jtMt4tIGvdlVm9V1ftF5JXE16GpuPLULIuqXuvdf9a756+o6rKEedfSaH8LcToM+Kqq3u6lfwLoVNX/NvJr1rNlNYZq6pemqqupfoeF+VIdf4XqnwOvVNUHxe0ZvNI7vQToU6NjQlXvSXj5qrR4a9EZx9uN/WG8hyymTK1SWFmW4Fpvt8rI+8s7gZl1lHFbSFiYwjqAktQhkLjOkpbli0A/8I068q6lhn8LCXUo7it3fR0I3Bg4XoHVhNN4/IsdAWxXVX8J+Kxx7kVgtohk/B+ziByAu1BqnKrS4rprSX6YtcrUkERkJq47+Ls6kr2I+66XU418pgEOcE6zyxhx/aR12JQ689pI/wE4Ut13FTcj79H8FuLKuxsuPJ82gg8Hvm4cH4r7ugWrCabxGJi9EvfHuq93bPZu3u+dPx9A3DF9P8dt49oB9HiN5eaPLyrtnriu089GWaZGNRt4a51p7gde8BvkRSQD/BrYt0VlDLt+0jocdXk82F+Nu7L6s17Y5QnzrvV7GM1vIU6HAct8eIv7OtvdcFcRN+M82oRrWY21Gm18BF6N+y/od6q8OnD+fFzLbh1wdeDcv3nnfgN8H9dludo7txfwey/PPwGneOFvAp7CXRr+NOATjCyzf3Eg7d3AX4CzQ8ryAeB9Rtrza5UJuCR4ncC9nOSVVXEb12/wtptw3byKcsaVxbuHPxj3/56E9RZZ3zW+s1NC4kTVYVU91Kiz0LKE1MXncF9SdYOxPRh3ryG/h6b+FkLq5AuMdNh8DrjcOHcisMI4TgGDwM5h6e3W3pt9hYCVVYz84Tyq+oXxSG81drJzma2srKw82aEBVlbxunOc01uNkazLbGVlZeXJusxWVlZWniwQraysrDxZIFpZWVl5skC0srKy8mSBaGVlZeXJAtHKysrK0/8HM9YLBEysRkEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "im.blur_circ(10*eh.RADPERUAS).display(\n",
    "                 plotp=True,pcut=0.1,mcut=0,\n",
    "                 label_type='scale',scalecolor='k',\n",
    "                 cfun='gray_r',cbar_unit=['Tb'],\n",
    "                 vec_cfun='rainbow',nvec=20,vec_cbar_lims=(0,1),scale_ticks=True,\n",
    "                 has_title=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1faea21",
   "metadata": {},
   "source": [
    "## Generate a synthetic observation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "676fe7f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate and save an EHT observation\n",
    "# ampcal and phasecal determine if gain variations and phase errors are included\n",
    "tint_sec = 60 #integration time in seconds, \n",
    "tadv_sec = 600 #advance time between scans\n",
    "tstart_hr = 0 #GMST time of the start of the observation\n",
    "tstop_hr = 24 #GMST time of the end of the observation\n",
    "bw_hz = 4.e9 #bandwidth in Hz\n",
    "obs = im.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz,\n",
    "                 sgrscat=False, ampcal=True, phasecal=False, ttype='direct')\n",
    "obs.save_uvfits(outpath+'.uvfits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8695920f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What is the resolution of the observation? \n",
    "beamparams = obs.fit_beam() # fitted beam parameters (fwhm_maj, fwhm_min, theta) in radians\n",
    "res = obs.res() # nominal array resolution, 1/longest baseline\n",
    "print(\"Clean beam parameters: [%.1f uas,%.1f uas,%.1f deg]\"\n",
    "      %(beamparams[0]/eh.RADPERUAS, beamparams[1]/eh.RADPERUAS,beamparams[2]/eh.DEGREE))\n",
    "print(\"Nominal Resolution: %.1f uas\"%(res/eh.RADPERUAS))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e34357",
   "metadata": {},
   "source": [
    "## Perform Total Intensity Imaging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6796e81e",
   "metadata": {},
   "outputs": [],
   "source": [
    "npix = 64\n",
    "fov = 160*eh.RADPERUAS\n",
    "zbl = im.total_flux() # total flux\n",
    "prior_fwhm = 100*eh.RADPERUAS # Gaussian size in microarcssec\n",
    "emptyprior = eh.image.make_square(obs, npix, fov)\n",
    "flatprior = emptyprior.add_flat(zbl)\n",
    "gaussprior = emptyprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0))\n",
    "gaussprior.display();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0b6c03",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Image total intensity\n",
    "imgr  = eh.imager.Imager(obs, gaussprior, prior_im=gaussprior, flux=zbl,\n",
    "                         data_term = {'amp':1, 'cphase':0.5},\n",
    "                         reg_term={'tv':0.05, 'l1':0.1},\n",
    "                         norm_reg=True, # this is very important!\n",
    "                         epsilon_tv = 1.e-10,\n",
    "                         maxit=250, ttype='direct')\n",
    "imgr.make_image_I(grads=True,show_updates=False)\n",
    "\n",
    "# Blur the image with a circular beam and image again to help convergance\n",
    "out = imgr.out_last()\n",
    "imgr.init_next = out.blur_circ(res)\n",
    "imgr.make_image_I(show_updates=False)\n",
    "\n",
    "out = imgr.out_last()\n",
    "imgr.init_next = out.blur_circ(res)\n",
    "imgr.maxit_next=500\n",
    "imgr.make_image_I(show_updates=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fdb0660",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Investigate output\n",
    "out=imgr.out_last().threshold(0.01)\n",
    "out.display();\n",
    "out.blur_circ(0.5*res).display();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3b486cf",
   "metadata": {},
   "source": [
    "## Perform Polarization Imaging on Fixed Stokes I background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98ba7a3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Image polarization with the polarimetric ratio as main data product\n",
    "imgr.init_next = out.blur_circ(0.25*res)\n",
    "imgr.transform_next = 'mcv'\n",
    "imgr.dat_term_next = {'m':1}\n",
    "imgr.reg_term_next = {'hw':1,'ptv':1}\n",
    "imgr.make_image_P(show_updates=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd13a636",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display Result\n",
    "out=imgr.out_last()\n",
    "out.display(plotp=True,pcut=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "967c97f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# blur and image again with the polarimetric ratio as data\n",
    "imgr.init_next = out.blur_circ(0,.5*res) # you need the extra argument to blur_circ to blur polarization!\n",
    "imgr.transform_next = 'mcv'\n",
    "imgr.dat_term_next = {'m':5}\n",
    "imgr.reg_term_next = {'hw':1,'ptv':1}\n",
    "imgr.make_image_P(show_updates=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f222f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display Result\n",
    "out=imgr.out_last()\n",
    "out.display(plotp=True,pcut=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a68e4a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# self-calibrate the data to the last image\n",
    "out=imgr.out_last()\n",
    "obs_sc = eh.selfcal(obs,out,ttype='nfft',use_grad=True,processes=-1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27a79da3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# image again using the self-calibrated data\n",
    "# image both total intensity and polarization simulatneously\n",
    "imgr.init_next = out.blur_circ(0,res) # you need the extra argument to blur_circ to blur polarization!\n",
    "imgr.obs_next = obs_sc # replace the observationwith the self-calibrated version\n",
    "imgr.transform_next = 'mcv'\n",
    "imgr.dat_term_next = {'amp':1, 'cphase':0.5,'pvis':1}  # you need data terms for both total intensity and pol\n",
    "imgr.reg_term_next = {'tv':0.05,'hw':1,'ptv':1} # you need regulraizer terms for both as well\n",
    "imgr.make_image_IP(show_updates=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f49bad05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display Result\n",
    "imgr.out_last().display(plotp=True,pcut=0.1);\n",
    "im.display(plotp=True,pcut=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fdccef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the output\n",
    "out.save_fits(outpath+'_output.fits')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
